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ABSTRACT 


This  report  summarizes  work  conducted  during  the  last 
year  under  the  contract  Research  in  Seismology.  These 
investigations  fall  within  the  four  broad  topics:  (1) 
Theoretical  and  observational  studies  of  source  mechanisms 
of  earthquakes  and  underground  nuclear  explosions;  (2) 

Wave  propagation  in  heterogeneous  and  anisotropic  media; 

(3)  Regionalization  of  earth  structure;  (4)  Array  and 
data  processing  techniques  for  event  detection  and  location 
All  publications  and  theses  completed  during  the  contract 
year  are  listed. 
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1 .  SUMMARY 


In  this  annual  report  we  review  work  completed  under 
the  contract  Research  in  Seismology  during  the  year 
1  April  1973  through  31  March  1974.  Within  the  broader 
problem  of  discriminating  earthquakes  from  underground 
nuclear  explosions,  we  have  conducted  a  number  of  specific 
investigations  of  seismic  source  and  wave  propagation 
studies . 

The  problems  considered  can  be  grouped  into  four 
major  subject  areas: 

(1)  Source  mechanisms  of  earthquakes  and  explosions 

(2)  Seismic  wave  propagation 

(3)  Earth  structure 

(4)  Arrays  and  data  processing 

In  the  following  sections  we  include  abstractr  of  papers 
in  each  category  either  published  or  soon  to  be  published. 
Recently  completed  work  is  discussed  in  greater  detail. 

We  also  list  all  publications  and  theses  supported  under 
this  project  during  the  contract  year. 

A  major  fraction  of  our  effort  was  devoted  to  continuing 
our  theoretical  and  observational  studies  of  the  source 
mechanisms  of  earthquakes  and  underground  nuclear  explosions. 
Work  by  Abe  on  two  earthquakes  in  Japan  and  by  Aki  and  co¬ 
workers  on  NTS  explosions  demonstrated  the  powerful  insight 
possible  by  combining  careful  near-field  and  far-field 
observations  of  seismic  sources.  A  very  important  result 
of  the  study  of  Aki  and  others  was  that  the  reduced 
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Jisplacement  pot-enti  for  nuclear  explosions  is  closer  to 
an  impulse  than  to  a  step  function,  providing  an  explana¬ 
tion  for  the  success  of  the  Mg"  m^  discrimination  criterion 
for  events  of  m,  >  4 . 

D  — 

Other  work  on  seismic  source  phenomena  included  Brown's 
study  of  premonitory  changes  in  seismic  wave  velocities  in 
the  source  region  of  strike-slip  earthquakes  in  Japan, 
continuation  by  Singh  and  Rosenman  of  their  calculations 
of  displacements  and  stresses  associated  with  various 
sources  in  a  viscoelastic  earth,  and  Richter  and  Johnson's 
novel  explanation  for  very  deep  earthquakes  as  associated 
with  standing  waves  on  a  mineralogical/chemical  transition 
boundary  at  650  km  depth  in  the  mantle. 

Several  wave  propagation  studies  were  completed.  Wang 

solved  the  problem  of  wave  propagation  in  weakly  anisotropic 

media  for  certain  geometries  of  wave  paths  and  symmetry 

directions,  with  direct  application  to  Rayleigh  and  P  waves 

n 

in  the  oceanic  upper  mantle.  Reasenberg  and  Aki  measured 
changes  in  the  in  situ  seismic  velocity  in  the  shallow  crust 
due  to  tidal  stresses.  Kuster  and  Toksflz  derived  the  elas¬ 
tic  moduli  of  a  two-phase  medium  using  a  wave-scattering 
approach,  and  applied  their  model  to  rocks  with  dry  and 
water-saturated  cracks. 

Much  of  our  attention  was  devoted  to  regionalization 
of  earth  structure,  particularly  to  the  study  of  areas 
within  and  adjacent  to  Asia.  In  his  M.S.  thesis.  Chapman 
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measured  the  dispersion  characteristics  of  Rayleigh  waves 
in  the  Arctic  Ocean  region  and  inferred  the  crustal 
structure  beneath  the  Arctic  Ocean  basin,  the  spreading 
Nansen  ridge  and  the  continental  shelves.  He  also  inter¬ 
preted  earthquakes  in  northeast  Asia  as  due  to  interaction 
between  the  Furasian  and  American  plates.  Solomon  and 
Butler  outlined  a  procedure  for  locating  seismically 
inactive  slabs  of  subducted  lithosphere  by  use  of  travel¬ 
time  residuals;  the  technique  may  be  useful  in  looking 
for  such  an  ancient  slab  propopcd  to  lie  beneath  the 
Urals.  Solomon  and  Julian  studied  the  effects  on  body- 
wave  propagation  paths  of  the  laterally  heterogeneous 
velocity  structure  beneath  active  spreading  centers.  Such 
effects  can  explain  the  observation  that  normal-faulting 
earthquakes  on  ridge  crests  appear  from  first-motion  dia¬ 
grams  to  be  accompanied  by  a  superimposed  'explosion' 
source. 

In  the  area  of  arrays  and  data  processing,  Shlien 
and  Toksfiz  extended  their  automatic  phase  identification 
algorithms  for  array  data  to  make  use  of  two  arrays 
running  concurrently.  The  techniques  were  tested  on 
LASA  and  NORSAR  records.  Ambuter  and  Solomon  built  a 
seismic  recording  package  using  an  event-recording  scheme 
based  on  a  continuously  updated  semiconductor  memory;  the 
scheme  results  in  far  less  storage  capacity  and  easier  later 
interpretation  than  do  conventional  continuous  recording 


methods . 


2.  SOURCE  MECHANISMS  OF  EARTHQUAKES  AND  EXPLOSIONS 

2<1  Seismic  Source  Function  for  a  1  Underground  Nuclear 

Explosion  by  Keiiti  Aki,  Michel  Bouchon,  and 

Paul  Reasenberg  (abstract) 

The  reduced  displacement  potential  obtained  from  close-in 
observation  of  seismic  displacement  during  an  underground 
explosion  usually  takes  the  form  of  a  step  function  with  or 
without  a  small  overshoot.  Theoretical  prediction  by  shock- 
wave  calculation  appears  to  agree  with  the  close-in  data. 

The  step-function  source  has  also  been  supported  by  the  ob¬ 
servations  on  Rayleigh  waves  at  periods  longer  than  10  sec. 

We  found,  however,  some  inconsistency  between  the  published 
data  on  residual  potentials  obtained  from  close-in  data  and 
those  on  seismic  moments  obtained  from  long-period  Rayleigh 
waves.  It  appears  that  only  about  1/3  of  the  residual 
potential  is  transmitted  to  the  far-field  at  long  periods. 

This  discrepancy  is,  however,  consistent  with  several  obser¬ 
vations  made  on  teleseismic  signals  suggesting  an  impulse 
rather  than  a  step  as  the  primary  form  of  the  potential 
function.  New  observations  of  the  two  NTS  events  at  dis¬ 
tances  2.6  to  7.8  km  using  wide  dynamic  range,  wide-band 
accelerometers,  combined  with  data  from  the  far-field, 
support  a  large  overshoot  4  to  5  times  the  residual  value. 

This  result  accounts  for  the  efficiency  of  the  Mg- 
discriminant  between  earthquakes  and  explosions  with  m. 
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around  4  and  greater.  The  compaction  of  the  source  volume 
by  spalling  was  suggested  as  a  possible  mechanism  for  the 
large  overshoot. 


2.2  Precursory  Changes  in  VP/VS  Before  Strike-slip  Events 
by  Raymon  Brown  (abstract) 

P  and  S-wave  travel  time  data  collected  in  southwestern 
Japan  by  means  of  a  five  station  microearthquake  array  are 
used  to  construct  Wadati  plots.  The  slopes  of  these  plots 
(VP/VS-1)  are  then  examined  as  a  function  of  time.  The 
largest  deviations  of  VP/VS  are  of  the  order  of  5% .  The  mag 
nitude  of  any  premonitory  changes  in  VP/VS  observed  here 
will  therefore  be  less  than  any  previously  reported.  Two 
of  the  largest  events  recorded  by  the  array  (M  *  5.0)  for 
which  there  is  enough  data  show  a  premonitory  increase  above 
the  mean  and  then  a  decrease  below  the  mean  and  finally  a 
rise  to  the  mean  shortly  before  the  event.  The  nature  of 
the  changes  before  both  events  agree  in  magnitude  and  time. 
The  time  from  the  minimum  VP/VS  to  the  time  of  the  events 
scales  in  agreement  with  previous  investigations.  The 
region  of  Japan  studied  in  this  paper  is  one  characterized 
by  shallow,  predominantly  strike-slip  events. 


2 • ^  Stress  Relaxation  in  a  Semi-Infinite  Viscoelastic 
Earth  Model  by  Martin  Rosenman  a  id  Sarva  Jit  Singh 
(abstract) 

Expressions  for  quasi-static  surface  stresses  resulting 
from  a  finite,  rectangular,  vertical,  strike-slip  fault  in 
a  Maxwellian  viscoelastic  half-space  are  derived.  Contour 
maps  are  obtained  in  some  representative  cases.  It  is  found 
that  all  nonvanishing  stress  components  at  the  free  surface 
die  exponentially  with  time.  This  is  in  contrast  to  the 
behavior  of  the  displacements  and  strains  which,  in  general, 
do  not  vanish  for  large  times. 

2*4  Stability  of  a  Chemically  Layered  Mantle  by  Frank  M. 
Richter  and  Carl  E.  Johnson  (abstract) 

A  system  of  two  immiscible  fluid  layers  and  heated 
from  below  is  analyzed  to  determine  its  stability  to  convec¬ 
tive  motions.  The  model  is  applicable  to  the  mantle  if 
some  fraction  of  the  density  increase  with  depth  results 
from  changes  in  chemical  composition .  For  different  para¬ 
meter  ranges,  three  modes  of  instability  are  possible: 

(1)  convection  of  the  entire  depth  of  the  fluid,  (2)  separate 
convection  within  each  layer,  and  (3)  overstability  in  the 
form  of  standing  waves  on  the  interface  between  layers. 

Earth  models  that  include  an  increase  in  density  due  to 
iron  enrichment  below  650  km  suggest  that  mantle  convection 


can  exist  above  and  below  but  not  through  the  650-km  level. 
The  model  results  are  combined  with  a  prior  understanding 
of  mineralogical  phase  changes  to  suggest  a  possible 
mechanism  for  the  deepest  earthquakes. 

2 . 5  Seismic  Displacement  and  Ground  Motion  Near  a  Fault: 

The  Saitama  Earthquake  of  September  21 ,  1931  by 

Katsuyuki  Abe  (abstract) 

The  dislocation  parameters  of  the  Saitama  earthquake 
(M  =  7.0,  36 . 15#N ,  139 . 24#E)  of  September  21,  1931,  are 
determined  on  the  basis  of  the  first-motion  data,  the 
aftershock  area,  and  the  close-in  seismograms  obtained  by 
a  low-magnification  seismograph.  The  earthquake  represents 
a  left-lateral  strike-slip  faulting  on  a  plane,  dipping 
80°  towards  N  196°E,  with  the  size  of  20  km  (length)  x  10  km 
(width) .  The  strike  of  the  fault  plane  is  found  to  be 
almost  parallel  to  that  of  the  eastern  extension  of  the 
Median  Tectonic  Line.  A  synthetic  study  suggests  that  the 
rupture  grows  bilaterally  at  a  velocity  of  2.3  km/sec. 

The  rise  time  and  the  final  dislocation  of  the  linear  ramp 
dislocation  time  function  are  determined  as  2  sec  and 
100  cm,  respectively.  The  corresponding  particle  velocity 
of  the  fault  dislocation  is  estimated  to  be  50  cm/sec. 

The  near-source  ground  displacements  and  ground  motions 
inferred  from  the  above  seismic  fault  model  are  found  to 


be  remarkably  consistent  with  the  data  on  the  high  precision 
levelling  and  on  the  field  survey  on  the  direction  to  which 
artificial  structures  collapsed. 


2 . 6  Fault  Parameters  Determined  by  Near-  and  Far-Field  Data: 

The  Wakasa-Bay  Earthquake  of  March  26,  1963  by 

Katsuyuki  Abe 
Summary 

The  source  process  of  the  Wakasa-Bay  earthquake  (M  =  6.9, 
35.80#N,  135.73 ’E,  depth  4  km)  which  occurred  near  the  west 
coast  of  Honshu  island,  Japan,  on  March  26,  1963,  is  studied 
on  the  basis  of  the  seismological  data.  Dynamic  and  static 
parameters  of  the  faulting  are  determined  by  directly  com¬ 
paring  synthetic  seismograms  with  observed  seismograms 
recorded  at  seismic  near  anl  far  distances.  The  De  Hoop- 
Haskell  method  is  used  for  the  synthesis.  The  average 
dislocation  is  determined  to  be  60  cm.  The  overall  dis¬ 
location  velocity  is  estimated  to  be  30  cm/sec,  the  rise 
time  of  the  slip  dislocation  being  determined  as  2  sec. 

The  other  fault  parameters  determined,  with  supplementary 
da  a  on  the  P-wa/e  first  motion,  the  S-wave  polarization 
angle  and  the  aftershocks,  are:  source  geometry,  dip 
direction  N  144°E,  dip  angle  68°,  slip  angle  22'  (right- 
lateral  strike-slip  motion  with  some  dip-slip  component) ; 
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fault  dimension,  20  km  length  x  8  km  width;  rupture  velo- 

25 

city,  2.3  km/sec  (bilateral);  seismic  moment,  3.3  x  10 
dyne  cm;  stress  drop,  32  bars.  The  effective  stress 
available  to  accelerate  the  fault  motion  is  estimated  to 
be  about  40  bars.  The  approximate  agreement  between  the 
effective  stress  and  the  stress  drop  suggests  that  most 
of  the  effective  s cress  was  released  at  the  time  of  the 
earthquake . 


1.  Introducticn 


One  of  the  most  important  factors  in  the  dynamic 
process  of  faulting  is  a  dislocation  velocity.  This  dynamic 
parameter  is  directly  related  to  the  effective  stress 
available  to  accelerate  the  fault  motion.  Such  effective 
stress  can  be  estimated  from  the  dislocation  velocity 
(e.g.  Kanamori,  1972).  Since  the  effective  stress  primarily 
controls  the  near-source  acceleration  (Brune,  1970),  its 
estimate  may  provide  the  possible  maximum  ground  acceleration 
expected  of  not  only  a  certain  earthquake  but  a  future 
earthquake.  As  the  determination  of  such  dynamic  parameters 
is  so  far  very  scanty,  detailed  study  on  the  fault  dynamics 
for  individual  earthquake  is  important  not  only  for  under¬ 
standing  the  physics  of  earthquakes  but  also  for  providing 
the  potential  data  for  earthquake  engineering. 

From  this  point  of  view,  we  will  attempt  to 
determine,  for  the  Waxasa-Bay  earthquake  of  March  26,  1963, 
the  dislocation  velocity  as  well  as  various  static  parameters 
such  as  fault  geometry,  fault  dimension  and  dislocation. 

The  near-field  seismograms  as  well  as  the  far-field  seismo¬ 
grams  are  available  for  this  earthquake.  Most  of  the  fault 
parameters  will  be  determined  by  means  of  a  direct 
comparison  between  synthetic  seismograms  and  observed 
seismograms . 
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2 •  The  Wakasa-Bay  earthquake  of  1963 
2.1 .  Redetermination  of  the  hypoconter 

The  Wakasa-Bay  earthquake  (M  =  6.9)  occurred  on 
March  26,  1963,  to  the  north  of  Kyoto,  near  the  west  coast 
of  Honshu  island,  Japan  (see  fig.  5).  The  hypocenters  of 
the  mam  shock  and  the  four  large  aftershocks  (m  =  4.8~5.2) 
were  redetermined  on  the  basis  of  the  P  times  as  reported 
in  the  Seismological  Bulletin  of  the  Japan  Meteorological 
Agency  ( JMA)  and  the  International  Seismological  Center 
(ISC)  for  the  year  1963.  The  Jeff reys-Bul len  travel-time 
table  was  used.  We  made  the  relocation  of  the  main  shock 
for  two  cases:  (1)  all  the  stations  with  epicentral 
distance  A  $  90°  were  used;  (2)  stations  with  A  £  11° 
were  used.  The  st xtions  used  in  the  case  (2)  belong  to  JMA. 
We  used  exclusively  tke  JMA  stations  for  the  four  aftershocks, 
since  the  world-wide  data  are  very  scanty.  in  all  cases 
the  threshold  of  the  O-C  residual  (observed  minus  computed 
P  time)  was  set  at  4  sec:  as  a  result  for  the  main  shock, 
throe  data  were  discarded  in  the  first  case,  and  five  data 
were  discarded  in  the  second  case.  The  number  of  stations 
included  in  each  set  is  shown  in  the  inset  of  fig.  i.  Also 
shown  is  the  number  of  stations  located  in  NE,  SE,  SW  and 
NW  quadrants  around  the  epicenter. 

The  root -mean-square  (RMS)  of  the  0-C  residuals 
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for  various  restrained  depths  are  shown  in  fig.  1  ,  where 
the  stability  of  the  depth  determination  is  shown.  The  RMS 
of  O-C  was  found  to  increase  as  the  depth  increased, 
favoring  a  shallow  focal  depth.  However,  it  is  well  known 
that  the  focal  depth  cannot  be  constrained  very  well  for  P 
tines  alone.  The  residuals  at  r.earby  stations  are  strongly 
affected  by  the  focal  depth.  This  feature  is  demonstrated 
in  fig.  2,  where  the  O-C  residuals  at  nearby  stations  with 
^  £  1° are  given  for  five  focal  depths.  It  is  evident  that 
the  O-C  residuals  increase  systematically  and  sharply  as 
the  depth  increased.  These  observations  lead  to  a  conclusion 
that  the  foci  of  the  main  shock  and  the  aftershocks  are 
very  shallow,  not  deeper  than  10  km.  In  the  following,  we 
will  use,  for  the  main  shock,  the  hypocenter  parameters 
determined  for  the  second  case:  origin  time  21**  34m  38.5s; 
latitude  35.80°N;  longitude  135.76°E;  depth  4  km.  Fig.  3 
shows  the  epicenters  of  the  main  shock  and  the  aftershocks 
which  occurred  within  twenty-four  hours  after  the  main 
shock.  The  epicenter  data  exc »pt  for  the  shocks  treated 
above  are  taken  from  the  bulletin  of  JMA. 

2.2.  Source  geometry 

The  source  geometry  of  the  main  shock  was  determined 
by  Ichikawa  (1971)  who  used  the  P-wave  first-motion  data 
reported  in  various  seisroological  bulletins,  and 


-13- 


tentativcly  by  Stauder  and  Bollinger  (1965)  who  used  mainly 
the  S-wave  polarization  angle  data.  As  the  hypocenter  is 
redetermined  here,  we  anew  determine  the  source  geometry. 

The  first-motion  data  of  P  waves  are  obtained  from 
the  long-period  seismograms  from  the  World  Wide  Standard 
Station  Network  (WSSN).  We  supplement  them  by  the  data 
reported  in  the  Seismological  Bulletin  of  JMA  and  the  data 
given  by  Stauder  and  Bollinger  (1965).  These  data  are 
plotted  on  the  Wulff  grid  as  shown  in  fig.  4.  The  lower 
hemisphere  of  the  focal  sphere  was  projected.  Though  two 
noc!?l  planes  can  be  defined  by  the  compression  and  dilatation 
fields,  we  introduce  the  S-wave  data  to  find  a  more  precise 
source-geometry.  The  15  data  of  polarization  angle  are 
taken  from  Stauder  and  Bollinger  (1965),  in  which  the 
polarization  angles  were  determined  from  particle  motion 
diagrams  of  S  phase  recorded  on  WWSSN  seismograms.  We  used 
Hirasawa ' s  (1966)  method.  The  focal  depth  of  the  earthquake 
is  placed  at  4  km.  The  two  nodal  planes  thus  determined 
are  remarkably  consistent  with  the  P  wave  data  as  shown  in 
fig.  4  by  solid  lines.  The  dip  direction,  the  dip  angle 
and  the  slip  angle  are:  <p  =  144°,  S  =  68°,  \  =  22°  for 
one  plane  and  f.  243%  £«  69°,  1  =  156°  fir  another 
plane.  The  standard  deviation  of  the  polarization  angles 
is  23*.  Although  some  inconsistent  data  exist  in  the 
P-wave  first  motions,  the  present  solution  is  believed  to 


-14- 


be  very  good  as  far  as  both  the  first  notion  and  the 
polarization  angle  data  are  concerned. 

Fxon  the  mechanism  diagram  alone,  it  is  not 
possible  to  select  the  actual  fault  plane  out  of  the 
orthogonal  nodal  planes.  The  spatial  distribution  of  the 
aftershock  provides  the  key  to  the  selection.  The  after¬ 
shocks  are  distributed  at  relatively  shallow  depths  and 
over  the  area  elongated  in  a  northeast-southwest  direction 
(fig.  3).  The  aftershocks  are  more  densely  distributed 
to  the  south-east  of  the  main  shock  epicenter  than  to  the 
north-west,  suggesting  the  fault  dipping  to  the  south-east. 

It  follows  that  the  actual  dislocation  took  place  over  the 
plane  striking  in  N  54°E  and  with  a  dip  angle  of  68*SE. 

The  size  of  the  fault  plane  is  estimated  as  20  km  length 
x  8  km  width  from  the.  extent  of  the  aftershock  area.  The 
present  solution  suggests  that  the  earthquake  represents 
an  almost  right-lateral  strike-slip  faulting  with  some 
dip-slip  component. 

^ •  Interpretation  of  near-field  seismograms 
3.1.  DataL_jand_me^hod__of  ana  lysis 

Records  of  a  near-source  ground  motion  were  obtained 
at  Maizuru  observatory,  50  km  away  from  the  epicenter,  and 
Abuyama  observatory,  106  km  away  from  the  epicenter.  Fig.  5 
shows  the  locations  of  the  earthquake  and  the  observatories 
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The  instrument  operated  at  Maizuru  observatory  is  a  strong- 
motion  seismograph  with  a  low  magnification:  its  constants 
are  always  fixed  to  be  To  (free  period  of  pendulum)  =6.0 
sec,  £  (damping  ratio)  =  8,  and  V  (static  magnification) 

=  1.0.  The  instrument  operated  at  Abuyama  observatory  is 
a  low-magnification  and  long-period  seismograph:  the 
seismogram  includes  a  calibration  pulse  from  which  the 
instrumental  constants  can  be  determined  as  To  =  28  sec, 

5.  =  2,3  and  ^  ~  1*1*  The  seismograms  are  shown  in  fig.  6. 

In  order  to  obtain  fault  parameters,  we  computed 
synthetic  seismograms  for  various  fault  models  and  compared 
them  with  the  observed  seismograms.  In  computing  the 
dynamic  near-field  displacement  the  integial  expressions 
given  by  Hat '  ell  (1969)  were  numerically  double-integrated 
over  a  fault  plane  (see  also,  De  Hoop,  1958).  It  is  assumed 
that  the  dislocation  takes  place  simultaneously  over  the 
fault  width  w  and  propagates  at  a  constant  velocity  v  along 
the  fault  length.  The  terapx ral  variation  of  the  dislo¬ 
cation  is  given  as  the  form  of  a  ramp  function  of  rise  time 
X.  ,  that  is 


f  0  t  <  0 

c(t)  =  i  t  /  t  o<t<r 
^  1  t  >  z 


(1) 
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Since  the  method  applies  only  to  an  infinite  homogeneous 
medium,  the  comparison  should  be  made  only  for  the  first 
several  cycles  of  body  waves.  A  free- surface  effect  was 
included  simply  by  doubling  the  amplitude  calculated  for 
an  infinite  medium.  Velocities  of  P  and  S  waves  are  placed 
at  6.0  km/sec  and  3.5  km/sec,  respect ively . 

3.2  Interpretation 

In  subsequent  calculations,  we  assume  several 
fault  parameters  as  follows:  144°,  6"  =  68°,  L  =  20  km, 

w  =  8  km.  From  fig.  3  the  rupture  is  considered  to  start  at 
the  center  of  the  fault,  and  to  propagate  bilaterally  along 
the  strike  of  the  fault.  With  these  constraint  we  can 
finally  control  the  rupture  velocity,  the  rise  time,  and 
the  final  dislocation'. 

From  the  tentative  examination  of  radiation 
patterns,  Maizuru  is  found  to  be  located  near  the  loop 
direction  for  S  waves,  and  near  the  node  direction  for  P 
waves.  On  the  o‘.her  hand,  Abuyama  is  found  to  be  located 
nearest  to  the  node  direction  for  S  waves:  the  predominant 
wave  following  P  waves  may  be  considered  to  be  surface 
waves,  probably  Rayleigh  waves,  in  view  of  the  late  arrival 
time.  On  these  grounds  we  use  only  the  S  wave  portion  for 
the  Maizuru  record,  and  the  P  wave  portion  for  the  Abuyama 
record,  for  a  reliable  estimate  of  fault  parameters. 
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Various  combinations  of  rise  times  (1,  2,  3,  4,  5 
«c)  and  rupture  velocities  (l.o,  l.S,  2.3,  3.0,  3.4  km/soc) 
are  assumed  in  an  attempt  to  deduce  the  most  probable  values. 
Synthetic  seismograms  are  given  in  fig.  7,  for  several 
combinations  of  rise  time  and  rupture  velocity.  A  rise  time 
and  a  rupture  velocity  affect  mainly  the  initial  slope 
and  duration  time,  respectively.  Comparison  with  the 
observed  data  suggests  that  combinations  of  ,  rise  time  of 
2  sec  and  a  rupture  velocity  near  2.3  km/sec  provide  an 
overall  agreement.  The  amount  of  the  dislocation  is  inde¬ 
pendent  of  the  nature  of  the  wave  form,  and  it  simply  scales 
the  amplitude.  The  observed  seismograms  are  -capered,  in 
amplitude,  with  the  synthetic  seismograms  caicu,ated  for 
t  =  2  sec  and  v  *  2.3  km/sec.  We  obtain  the  dislocation 
of  62  cm  and  68  cm  fot  the  Maizuru  and  Abuyame  records, 
respectively;  the  average  is  65  cm. 


Although  the  dislocation  simply  scales  the  amplitude, 
the  amplitude  is  very  sensitive  to  both  the  rupture  velocity 
and  the  rise  time.  For  example,  the  amplitude  tends  to 
increase  with  the  rupture  velocity  increased,  as  shown  in 
fig.  7.  if  the  dislocation  can  be  independently  determined, 
the  present  estimate  of  v  and  Z  may  be  fairly  confirmed. 

From  this  reason  we  try  to  analyze  far-field  seismograms  to 
determine  the  dislocation  independently. 
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4 .  Interpretation  of  far-field  seismograms 
4.1.  Data  and  method  of  analysis 

The  copies  of  the  long-period  seismograms  from  the 
WWSSN  stations  were  collected.  We  used  horizontal  components 
for  7  S-waves ,  and  vertical  components  for  12  P-waves . 

Fig.  8  show*  the  theoretical  radiation  pattern  of  P,  SH  and 
SV  waves,  and  the  azimuthal  distribution  of  the  stations 
used.  It  is  to  be  noted  in  fig.  8  that  the  SH  motion  excels 
the  SV  motion  in  amplitude  at  a  few  stations  located  in 
both  the  north-east  and  north-west  directions.  With  this 
geographical  filter  we  exclusively  used  the  S  waves  pre¬ 
dominated  by  SH  wavo«,  in  order  to  avoid  ar.y  possible 
contamination  which  may  result  from  reflection  and  refraction. 
The  stations  and  other  pertinent  data  used  here  are  listed 

i 

in  Table  1 . 

For  the  interpretation  of  the  seismograms ,  synthetic 
seismograms  are  computed  for  various  propagating  fault 
models.  For  the  computation,  we  employed  the  Haskell's 
(1964)  method  which  gives  the  far-field  displacement  in  an 
infinite  homogeneous  medium.  The  setup  of  modeling  a  fault 
source  is  described  in  the  preceding  section.  In  the 
framework  of  this  model,  the  far-field  displacement  can  be 
expressed  in  a  good  approximation,  as  the  form  of  a 
trapezoidal  pulse  by 
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Uc  *  C  (W  *<t;  X  ’  *-) 


(2) 


where 


C  = 


el 


££_ 


4-710  y  3 

->  c 


(3) 


f(t;  T  ,  t  )  =  G(t  -  —  )  G(t - —  -  t  ) 

T/c  vc  ^ 

tc  =  v(1  ’  IF  c<‘®*  >  (5) 


(4) 


In  these  expressions  U£  =  far-field  displacement  (suffix  c 


denotes  to  P,  SH,  SV  waves),  Mq  =  seismic  moment,  - 

radiation  pattern,  j5  density,  v  =  wave  velocity,  r  = 

C  O 


dlst 3incA  < 


♦  />•  •  _  _  _  _  1  _  »  .  . . 

uwm4.wc  —  ttuvjit:  D«r  i  Vv i  t lie  lay 

path  and  rupture  direction.  The  displacement  source  time 

,  .  .  Equation  4  is  illustrated  in  fig.  9 . 

function  G(£)  is  defined  in  eq .  1.^  Equation  2  is  readily 

generalized  to  bilateral  rupture  by 

"c  *  <c/2)  Mo  [  f(li  Z  .  tb  /  t* 


f(‘i  T  .  t‘)  /  t* 


] 


(6) 


where 


=  2^  (1  ’  ^  c°s  *  > 
*2  5  » 


(7) 

(8) 


For  a  direct  comparison  between  synthetic  and 
observed  seismograms,  various  corrections  must  be  included. 
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Such  corrections  take  care  of  the  geometrical  spreading, 
the  instrumental  response,  the  amplitude  and  phase  distortions 
due  to  a  layered  crust  (Haskell,  1960,  1962),  and  the 
amplitude  and  phase  distortions  due  to  the  anelasticity  of 
the  Earth  (Futterman,  1962).  All  these  responses  were 
included  c.  frequency  domain.  The  instrumental  response  was 
calculated  from  the  Hagiwara’s  (1958)  formula;  the  parameters 
are  To  =  30  sec,  Tg  =  100  sec,  hp  =  1 ,  hg  .  t (  and  ^ 

(coupling  factor)  =  0.1.  llie  attenuation  factor  and  the 
crustal  structure  are  the  same  as  those  used  by  Fukao  (1970). 
The  S  waves  are  synthesized  from  both  SH  and  SV  waves. 


4-2«  Interpretation 

For  the  computation  of  synthetic  seismograms,  we 
use  eq.  6.  In  the  present  study  take-off  angles  of  the  ray 
paths  range  from  15“  to  22°  for  both  P  and  S  waves.  Thic 
small  range  of  take-off  angles  and  the  poor  azimuthal 
coverage  of  the  stations  are  not  good  enough  to  adequately 
determine  the  rupture  velocity.  Then  the  rupture  velocity 
is  constrained  at  2.3  km/sec.  In  these  cases,  the  rupture 
term,  (v/v^)  cos  0  ,  in  eq.  7  and  8  is  evaluated  to  be  less 
than  0.1  to  0.2.  This  term  is  very  small  as  compared  with 
unity,  so  that  we  may  neglect  them.  Then,  we  can  simply 
place  Ctc>  -  tc  in  eq.  6.  Since  the  wave  form  at  far 

field  is  not  very  sensitive  to  the  rise  time  because  of  the 
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anelast icity  of  the  medium,  we  reasonably  assume  T  =  2  sec. 

With  these  parameters  we  can  compute  synthetic  seismograms 

so  as  to  vary  <t  >  and  M  The  parameters  <t  >  and  M 
—  c  o 

affect  mainly  the  width  of  the  initial  pulse  and  the~ampli- 
tude,  respectively. 

A  direct  comparison  of  the  best  fit  synthetic 
seismograms  with  the  observed  seismograms  is  shown  in  fig.  10 
for  S  waves  and  in  fig.  11  for  P  waves.  In  the  present 
study  it  is  attempted  to  fit  only  the  first  half  cycle.  It 
can  be  seen  in  terms  of  <t^>  and  the  synthetic  seismograms 
explain  reasonably  the  observed  wave  forms  of  at  least 
the  first  half  cycle  of  both  S  and  p  waves.  The  discrepancy 
after  the  first  half  cycle  of  the  pulse  arrival  is  probably 
due  to  the  effect  of  the  free  surface  near  the  source. 
Actual*/  a  good  agreement  of  general  features  between  the 
observed  and  synthetic  wave  forms  has  been  obtained  at 
least  for  deep  earthquakes  (e.g.,  Fukao,  1970;  Mikumo,  1971). 

By  a  direct  comparison  between  the  synthetic  and 
observed  seismograms  for  the  pulse  width  and  the  amplitude 
of  the  first  half  cycle  of  the  S  and  P  waves,  we  obtain  the 
results  for  <t^>  and  (Table  1).  The  seismic  moment  is 
determined  as  2.6  x  1025  dyne-cm  for  S  waves  and  3.3  x  1025 
dyne-cra  for  P  waves;  the  average  value  is  3.0  x  1025  dyne-cm. 
According  to  the  slip  dislocation  theory  of  the  faulting 
(Aki,  1966),  the  slip  dislocation  <D>  averaged  over  the 
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fault  area  can  be  readily  estimated  from  <D>  =  M  /jjlS, 

where  yx  is  the  rigidity  and  S  is  the  fault  area.  Using 

the  values  Mp  =  3.0  x  10  **  dyne-cm,  JJ-  =  3.4  x  10^  dyne/cra^, 
2 

S  =  20  x  8  km  ,  we  have  <D>  =  55  cm. 

We  can  also  estimate  the  fault  length  from  <t  >. 

We  obtain  L  =  17  +  8  km.  The  poor  accuracy  of  the  result 
is  partly  due  to  the  simplification  of  the  rupture  term. 

The  point  to  be  emphasized  here  is  that  the  fault  length 
estimated  from  the  aftershock  area  is  not  significantly 
different  from  that  estimated  from  the  seismic  waves. 

^ •  Dislocation  velocity  and  effective  stress 

We  obtained  <D>  =  65  cm  from  the  seismic  near-field 
data  and  <D>  =  55  cm  ^frcm  the  teleseismic  data,  independently. 
We  consider  that  the  agreement  is  reasonably  gr»od.  As 
described  before,  this  close  agreement  gives  us  confidence 
in  the  previously  obtained  results,  in  particular,  X  *  2 
sec  and  v  =  2.3  km/sec.  The  amount  of  dislocation,  60  cm, 
is  used  in  the  following  discussion. 

The  dislocation  velocity  (the  velocity  of  one  side 
of  the  fault  with  respect  to  the  other)  can  be  determined 
from  the  dis location-time  history  during  the  earthquake. 

Since  the  rise  time  which  represents  the  time  required  for 
the  completion  of  the  dislocation  at  a  point  on  the  fault 
has  been  obtained,  we  can  determine  the  overall  dislocation 
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velocity  dividing  th"1  final  dislocation  by  the  rise  time, 
with  <D>  =  60  cm  and  Z  -  2  sec,  the  result  is 

<D>  =  <D>  /  Z  =30  cm/sec  (9) 


Then,  the  overall  particle  velocity  of  the  faulting  (one-half 
the  velocity  of  one  side  of  the  fault  with  respect  to  the 
other)  is  estimated  to  be  15  cm/sec. 

The  dislocation  velocity  is  directly  related  to 
the  effective  stress  which  accelerates  the  fault  motion. 
Kanamori  (L972)  has  related  the  overall  dislocation  velocity 
with  the  effective  str  ss,  summarizing  the  presently 


available  dvnami c-di O nrat \ on 
According  to  his  summary,  the 


models  of  the  faulting, 
most  practical  expression  for 


a  bilateral  fault  of  finite  length  and  a  ^i-.jite  rupture 

I 

velocity  is  given  by 


(1  ♦ 


A  A  <D> 

U  1 3  2 


(10) 


where  rfeff  is  the  effective  stress,  the  shear  velocity, 

v  the  rupture  velocity,  jl  the  rigidity  and  <D>  the 
overall  dislocation  velocity.  Using  =  3.5  km/sec,  = 

2.3  km/sec,  fJL  =  3.4  x  1011  dyne/cm*  <D>  =  30  cm/sec,  we 
obtain  the  effective  stress  of  37  bar. 

Using  the  slip  dislocation  theory  of  faulting 
(Knopoff,  1958),  we  can  estimate  the  static  fault  parameters 
such  as  the  stress  drop  £  ,  the  strain  drop  ,  and 
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the  released  strain  energy  W.  The  formulas  used  here  are: 

„  v  4  fi<D> 

AO  ^  —  £- -  ,  At  = 

71  w  1 

W  =  -2-  jul  < D>2  L 

71  r 

The  results  are:  A  &  -  32  bar;  =  0.94  x  10  4;  W  = 

21 

1.6  x  10  erg.  The  stress  drop  is  about  the  same  order  as 
that  for  large  shallow  earthquakes  which  occurred  beneath 
the  trench  along  the  Pacific  coast  (e.g.  Abe,  1972).  The 
stress  drop  refers  to  the  difference  in  stress  before  and 
after  the  formation  of  the  fault,  but  not  to  the  prevailing 
tectonic  stre*;*;  Hnwovor,  it  is  to  bo  noted  that  the  stress 
drop  is  about  the  same  as  the  effective  stress.  This 
approximate  equality  suggests  that  most  of  the  effective 

I 

stress  which  had  been  acting  on  the  fault  prior  to  the 
earthquake  was  released  at  the  time  of  the  earthquake. 

6.  Discussion  and  conclusion 

Some  major  earthquakes  occurred  in  the  neighboring 
regions  of  the  Wakasa-Bay  earthquake  (M  =  6.9)  of  1963. 

Major  destructive  shocks  are  the  Tottori  earthquake  (M  = 

7.4)  of  1943  and  the  Fukui  earthquake  (M  =  7.3)  of  1948. 
Kanamori  (1972,  1973)  obtained,  from  the  analysis  of 
close-in  seismograms,  the  overall  dislocation  velocity  of 
83  cm/sec  for  the  Tottori  earthquake  and  lOOcra/sec  for  the 
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Fukui  earthquake.  The  two  earthquakes  represent  a  predomi¬ 
nant  strike-slip  faulting.  We  note  that  these  estimates  are 
about  three  times  as  large  as  the  dislocation  velocity  of 
the  Wakasa-Bay  earthquake.  In  other  words,  it  is  such  a 
difference  of  the  overall  dislocation  velocity  that  suggests 
a  significant  difference  in  the  stress  level  in  the 
respective  epicentral  areas. 

We  further  note  that  the  earthquakes  cited  above 
represent  almost  the  complete  release  of  the  effective 
stress  at  the  time  of  the  earthquake.  This  feature  is  what 
might  be  expected  of  a  purely  elastic  rebound  in  which  the 

S  t  Or  Pfl  QTKOQQ  1  C  f  ol  03  CoH  nc  tint  nnn  A(al  •  ^  ^ 
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brittle  fracture.  In  such  a  case,  the  stress  drop  is 
indicative  of  the  effective  tectonic  stress  and  then  the 

I 

dislocation  velocity  can  be  eventually  interpreted  in 
terms  of  the  stress  drop.  By  replacing  eq .  1 ^  by 

<6>  ^  -£-  )-l*cr  (i2) 

H-  v 

we  can  estimate  the  overall  dislocation  velocity.  In  the 
limit  of  y  -  p  ,  as  a  rough  approximation,  eq .  12  gives 

<*>>  **  (  ~~  )  (13) 

r 

This  relation  may  be  useful  for  predicting  the  overall 
dislocation  velocity  expected  for  a  future  large  earthquake, 
for  example,  if  stress  drops  are  determined  for  small 
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earthquakcs  occurring  at  the  expected  localities. 

Oi  the  basis  of  the  seismic  near-field  seismograms, 
the  far-field  seismograms  and  the  aftershocks,  the  source 
process  of  the  Wakasa-Bay  earthquake  of  March  26,  1963,  are 
studied.  Synthetic  seismograms  were  computed  for  various 
fault  models  and  were  compared  with  the  observed  seismograms 
to  obtain  appropriate  fault  parameters.  The  results  can  be 
summarized  as  follows. 

Hypocenter:  origin  time  21 h  34m  38.5s,  latitude 
3 5 . 80°  N ,  longitude  135.76°E,  depth  4  km, 

Fault  geometry:  dip  direction  N  144°E,  dip  angle  68°, 

M  1  «'  «  ,  i  o 

ouyi^?  4 c*.  t 

Fault  motion:  right- lateral  strike-slip  motion 
with  some  reverse  dip-slip  component, 

Fault  dimension:  length  20  km,  down-dip  width  8  km, 

Rupture  velocity:  2.3  km/sec  (bilateral), 

25 

Seismic  moment:  3.3  x  10  dyne-cm, 

Average  dislocation:  60  cm, 

Stress  drop:  32  bar, 

-4 

Strain  drop:  0.94  x  10  , 

Rise  time:  2  sec, 

Overall  dislocation  velocity:  30  cm/sec, 

Effective  stress:  40  bar. 

It  is  concluded  that  the  earthquake  represents  a  nearly 
complete  release  of  the  effective  stress.  For  such  a  case, 


the  overall  dislocation  velocity  is  approximately  given  by 
the  relation,  <D>  %  (  6  /  v  )  A  o,  where  B  is  the  shear 

wave  velocity,  y  the  rigidity,  and  A  o  the  stress  drop. 
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Fig.  1  The  root-mean-square  (RMS)  of  the  observed  minus 

computed  P  times  (O-C  residuals)  versus  focal  depth. 

The  total  number  of  stations  in  the  NE,  SE,  SW  and 
W  quadrants  is  shown  in  the  inset. 

Fig.  2  O-C  residuals  at  nearby  stations  for  5  focal  depths. 
Approximate  distribution  of  7  stations  (A<  1°)  is 
shown  on  the  bottom. 

Fig.  3  Distribution  of  the  aftershocks  within  24  hours  after 
the  main  shock. 

Fig.  4  The  P-wave  firrt  motion  and  the  S-wave  polarization 
angle  obtained  for  the  Wakasa-Bay  earthquake  o.r 
March  26,  1963.  The  lower  half  of  the  focal  sphere 
is  projected  on  the  Wulff  grid.  f  is  the  dip  direction 
measured  clockwise  from  the  north  and  S  is  the  dip 
angle.  The  solid  curves  show  the  nodal  plane  solution 
determined  by  the  S  wave  data. 

Fig.  5  The  fault  geometry  of  the  Wakasa-Bay  earthquake  of 
1963.  The  locations  of  the  Maizuru  and  Abuyama 
seismological  observatories  are  shown. 

Fig.  6  Seismograms  of  the  Wakasa-Bay  earthquake  of  1963 
recorded  at  Maizuru  and  Abuyama.  To  is  the  free 
period  of  pendulum,  £  the  damping  ratio  and  V  the 
static  magnification. 

Fig.  7  Observed  seismograms  and  synthetic  seismograms 
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calculatcd  for  various  values  of  rise  time  T  and 
rupture?  velocity  \/.  The  dislocation  <D>  is  held 
constant . 

Fig.  8  Radiation  patterns  of  P,  SH  and  SV  waves  for  the  take¬ 
off  angle  of  20*.  The  azimuthal  distribution  of 
the  stations  used  is  shown. 

Fig.  9  General  form  of  the  source  time  function  f(t;  %  »t  ) 

in  the  Haskell  model.  The  cases  for  T<  t  and 

c 

^  ^  are  shown.  Both  areas  are  equal  to  t  . 

Fig.  10  Comparison  of  the  observed  far-field  seismograms  of 
S  waves  with  the  best  fit  synthetic  seismograms. 
Amplitude  is  normalized. 

Fig.  11  Comparison  of  the  observed  far-field  seismograms  of 
P  waves  with  the  best  fit  synthetic  seismograms. 
Amplitude  is  normalized. 


Table  1.  Station  data,  pulse  width  and  seismic  moment 


I 


h  N  o  o  O  o  o  m  cn  cn 

•  •••i  ■  •  . . 

pi  ^  n  V  cn  cn  cn  cm  cn  cn  cn 


cn  cn  cn  cn  i^f^r  i  cm  cm  cn  cm  m  cm 


O' 

H 

m 

H 

<r 

• 

•  1 

i  •  i 

1 

1 

• 

(V 

CM 

CM 

CM 

i  r>-  i  *r  t*- 


i  >  i  i 


rf  in  m 


iD'OOcoo'CDifi'Ot^coina'OO 
f— 4  H  PM  r~t  H  W  »—4  *— 4  H  H  W  PM 


oo 

m 

cn 

O' 

CM 

cn 

O' 

CM 

cn 

in 

H 

vO 

o 

O' 

CM 

o 

O' 

CM 

CO 

*t 

cn 

-o 

m 

O 

CM 

Tf 

o 

t" 

cn 

CM 

CM 

*t 

CM 

o 

r* 

H 

m 

'T 

m 

H 

cn 

cn 

cn 

CM 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

cn 

M. 

•  •k 

u 

-H 

4-* 

0) 

C 

a 

•  M 

0) 

V 1 

£ 

B 

01 

«-* 

0 

u 

3 

B 

B 

•k 

•H 

U 

c n 

N 

•H 

0) 

<« 

B 

> 

in 

<0 

•H 

3 

0)  o  c 

o  Z 


CM 

o 

CM 

O' 

CM 

CM 

Q 

■o 

cn 

cn 

H 

CO 

cn 

CO 

CO 

t" 

o 

H 

CM 

CM 

o 

CM 

CO 

cn 

in 

00 

cn 

o 

o 

CM 

m 

-o 

CM 

CO 

m 

CO 

co 

cn 

CO 

CO 

CO 

\D 

CO 

CO 

CO 

o 

CO 

o 

H  4$ 

X3  3 


ci  ^ 

o 


u 

H  2  2 

:  ce 

3  CD 

3 

Cfl  O  2 

H  X 

(-4*5 

l  2 

V)  5 

width  of  S 


FUKUI 


3*«9-01) 

nanzivw 


CCG 


-45- 


3.  SEISMIC  WAVE  PKOPAGATION 

3.1  Elastic  Wave  Propagation  in  Homogeneous  Transversely 
Isotropic  Medium  With  Symmetry  Axis  Parallel  to  the 
Free  Surface  by  Chi-fung  Wang  (abstract) 

Phase  velocities  of  elastic  wave  propagation  in  a 
homogeneous  transversely  isotropic  medium  with  symmetry  axis 
parallel  to  the  free  surface  of  a  half  space  is  investigated. 
Approximate  solutions  of  the  problem  of  pha^e  velocities  of 
Rayleigh,  horizontally  propagating  P  and  SH  waves  is  obtained 
by  means  of  perturbation  method  on  the  assumption  that  the 
deviation  of  the  elastic  coefficients  from  isotropy  is  small. 
In  the  case  of  horizon  ally  propagating  SV  waves  an  exact 
solution  is  obtained.  The  vertical  lamination  model  approxi- 
mating  fracture  zones  and  the  Olivine  model  showing  large 
azimuthal  variation  of  P  and  Rayleigh  waves  needed  some 
modification  in  such  a  way  that  the  a  axis  of  Olivine  crystal 
will  be  distributed  diffusely.  Suitable  choice  of  weighting 
functions  averaging  the  orientation  of  an  a  axis  will  give 
close  agreement  between  both  observed  and  predicted  P  and 
Rayleigh  waves. 
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3.2  A  Precise,  Continuous  Measurement  of  Seismic  Velocity 
for  Monitoring  In  Situ  Stress  by  Paul  Reasenberg  and 
Keiiti  Aki  (abstract) 

An  air  gun  repeatedly  shot  in  a  water-filled  hole  every 
6  or  10  s  was  used  to  measure  the  in  situ  seismic  velocity 
over  a  distance  of  200  m  in  a  granite  quarry.  We  found  a 
P«ak-to^p«ak  half-percent  variation  of  velocity,  which 
correlates  wsll  with  the  variation  in  tidal  stress.  The 
inferred  high  stress  sensitivity  of  velocity  change  (0.2 
bar  *)  nay  be  attributed  to  the  presence  of  extremely  thin 

-5 

cracks  (aspect  ratio  10  )  by  which  the  rock  mass  near  the 

surface  is  broken  into  blocks .  Such  cracks  should  close  at 

* 

a  depth  of  about  15  m.  When  the  present  results  for  de¬ 
tecting  the  tidal  effect  on  velocity  are  extrapolated  for 
waves  penetrating  to  depths  of  a  few  kilometers,  the  velo¬ 
city  must  be  measured  with  a  precision  of  better  than  one 
part  in  105. 


3 . 3  Velocity  and  Attenuation  of  Seiimic  Waves  in  Two-Phase 
Media:  I.  Theoretical  Formulations  by  Guy  T.  Kuster 

and  M.  Nafi  Toksflz 
Summary 

The  propagation  of  seismic  waves  in  two-phase  media 
is  treated  theoretically  to  determine  the  elastic  moduli 
of  the  composite  medium  given  the  properties,  concentration 


and  shapes  of  the  inclusions  and  the  matrix  material. 
Assuming  long  wavelengths  the  problem  is  formulated  in  terms 
of  scattering  phenomena  similar  to  the  approach  of  Ament 
(1959).  The  displacement  fields,  expanded  in  series,  for 
waves  scattered  by  an  "effective"  composite  medium  and 
individual  inclusions  are  equated.  The  coefficients  of 
series  expansions  of  the  displacement  fields  provide  a 
relationship  between  the  elastic  moduli  of  the  effective 
medium  and  those  of  the  matrix  and  inclusions.  The  expres¬ 
sions  are  derived  both  for  solid  and  liquid  inclusions  in 
solid  matrix,  as  well  as  for  solid  suspensions  in  fluid 
matrix.  Both  spherical  and  oblate  spheroidal  inclusions 
are  considered. 

Some  numerical  calculations  are  carried  out  to  demon¬ 
strate  the  effects  of  fluid  inclusions  and  inclusion  shapes 
on  the  seismic  velocities  in  rocks.  It  is  found  that  in 
addition  to  concentration,  inclusion  shapes  and  properties 
are  important  parameters.  A  concentration  of  a  fraction  of 
one  percent  of  thin  (small  aspect  ratio)  inclusions  could 
effect  the  compressional  and  shear  velocities  by  more  than 
ten  percent.  For  both  sedimentary  and  igneous  rock  models, 
the  calculations  for  "dry"  (i.e.  air-saturated)  and  water- 
saturated  states  indicate  that  the  compressional  velocities 
change  significantly  while  the  shear  velocities  change 
much  less  upon  saturation  with  water . 
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LIST  OF  SYMBOLS 


A  amplitude  of  incident  plane  P  wave 

A ,  ,  fourth  order  tensor 

ijpq 

a  radius  of  spherical  inclusion 

B  coefficient  in  die  series  expansion  of  scattered 

n  P  waves 

j  volume  concentration  of  inclusions 


ijpq 

:ijpq 


ij 

:ki 

M  \ 

C* 


j 


n 


K 

K' 


elastic  tensor,  matrix 

elastic  tensor,  inclusion 

strain  field 

incident  strain  field 

Green's  function,  matrix 

spherical  Hankel  function,  first  kind 

spherical  Bessel  function 

matrix  bulk  modulus 

inclusion  bulk  modulus 


K*  effective  bulk  modulus 


1 


ij 


N 


direction  cosine 

number  of  inclusions  in  representative  sphere 


p  wave  number,  P  wave  in  matrix 

pi  wave  number,  P  wave  in  inclusion 

P  (cos 9)  Legendre  polynomial,  order  n 
n 

r  radius,  representative  sphere 


s 


wave  number,  S  wave  in  matrix 
wave  number,  S  wave  in  inclusion 


Tijpq 

uijpq 


u 


u* 

Au 

u 

u* 

V 


v 


vi 

v* 


X 


*0 


4th  order  tensor  depending  on  matrix  and  inclusion 
properties 

fourth  order  tensor 

total  displacement  vector 

incident  displacement  vector 

scattered  displacement  vector  for  s^h  inclusion 
scattered  displacement  vector  from  representative  sphere 
scattered  displacement  vector 
radial  displacement 

radial  displacement  scattered  from  jth  inclusion 

radial  displacement  scattered  from  representative  sphere 

volume  of  an  inclusion 

volume  of  representative  sphere 

volume  of  inclusion  in  representative  sphere 

displacement  vector  inside  an  inclusion 

transverse  displacement 

transverse  displacement  scattered  from  inclusion 

transverse  displacement  scattered  from  representative 
sphere 

point  in  the  matrix 

center  of  representative  sphere 

center  of  8th  inclusion 
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<x  aspect  ratio  of  oblate  spheroid 

a*  effective  P  wave  velocity 

a*  effective  "static"  P  wave  velocity 

6*  effective  S  wave  velocity 

X  matrix  Lame's  constant 

X'  inclusion  Lamp's  constant 

li  matrix  shear  modulus 

P'  inclusion  shear  modulus 

P*  effective  shear  modulus 

P  matrix  density 

P '  inclusion  density 

p*  effective  density 

Pj  effective  gravitational  density 

L  point  inside  an  inclusion 

£.  center  of  an  inclusion 

£°  center  of  representative  sphere 

n  viscosity 

u  angular  frequency 
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INTRODUCTION 

In  the  earth,  where  rocks  are  generally  saturated  or 
partly  saturated  with  fluids  the  seismic  waves  propagate 
through  two-phase  media.  For  theoretical  treatment,  a 
two-phase  medium  is  defined  as  an  aggregate  of  two  homo¬ 
geneous  phases  of  different  properties,  where  one  phase 
(the  matrix)  is  a  continuum  in  which  inclusions  of  the  other 
phase  are  randomly  embedded.  If  the  two-phase  medium  is 
quasi-homogeneous ,  one  could  define  an  homogeneous  medium 
(the  effective  medium)  which  is  equivalent  to  the  two-phase 
medium  on  a  macroscopic  scale.  In  this  paper  we  will 
derive  theoretical  expressions  for  the  effective  properties 

of  a  two-phase  medium  for  the  propagation  of  elastic  waves 

% 

whose  wavelengths  are  much  longer  than  the  size  of  an 
inclusion. 

Theoretical  treatment  of  the  propagation  of  elastic 
waves  in  two-phase  media  (the  dynamic  problem)  is  relatively 
scarce-  In  a  few  studies  on  this  subject  (Ament,  1959;  Mai 
and  Knopoff ,  1967) ,  the  results  were  restricted  to  the  case 
where  the  matrix  is  solid,  the  inclusions  are  spherical  and 
much  smaller  than  the  wavelengths  and  sufficiently  far  apart 
from  each  other  so  that  interactions  are  negligible.  Ament 
(1953)  studied  the  case  of  rigid  spheres  in  a  viscous  liquid 
matrix.  Biot  (1956a, b)  treated  the  problem  for  porous  rocks 
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for  both  high  and  low  frequency  limits.  His  elegant  formulations 
require  the  specification  of  a  number  of  parameters  before 
the  effective  medium  properties  can  be  computed.  Chekin  (1970) 
worked  on  the  problem  of  wave  propagation  in  rocks  with 
cracks,  but  his  formulations  were  for  zero-width  cracks. 

The  theoretical  elastic  behavior  of  two-phase  media  under 
static  loading  has  been  studied  in  detail  (see  Hashin,  1970, 
for  a  review) .  The  results  of  these  studies  have  been  applied 
to  seismic  problems  under  the  assumption  that  the  conditions 
prevailing  in  the  propagation  of  long  wavelength  waves  can 
be  approximated  by  those  under  static  loading  (Eshelby,  1957; 
Hashin,  1962;  Hashin  and  Shtrikman,  1963;  Wu,  1966;  Walsh, 

1969;  Solomon,  1972;  Dederichs  and  Zeller,  1973;  Korringa, 

1973;  Zeller  and  Dederichs,  1973).  The  static  approaches, 
although  not  exact  formulations  for  the  dynamic  problem  have 
lead  to  most  useful  results.  As  our  comparisons  will  show 
later,  scattering  phenomena  and  inertia  effects,  intrinsic 
features  o.*  4,*e  wave  propagation  generally  omitted  in  static 
models,  become  of  importance  only  in  limited  cases. 

In  this  paper  we  follow  the  approach  of  Ament  (1959) 
and  formulate  the  problem  in  terms  of  scattering  phenomena. 

We  repeat  his  derivation  for  spherical  inclusions  embedded 
in  a  solid  matrix  and  we  treat  the  additional  cases  of  a 
fluid  matrix  and  of  spheroidal  inclusions.  All  our  models 
involve  the  assumptions  that  th'i  wavelengths  are  much  longer 
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than  the  size  of  the  inclusions  and  that  multiple  scattering 
effects  are  negligible. 

THEORETICAL  FORMULATION  IN  TERMS  OF  SCATTERING 

Consider  N  inclusions  randomly  embedded  within  a  finite 
region  VQ  of  an  infinite  matri.:  (Figure  1)  .  Let  an  elastic 
wave  be  incident  from  infinity.  We  may  write  tne  displacement 
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u(x)  observed  at  a  point  x  outside  as 


(1) 


Q 

where  u  (x  >x^)  is  the  displacement  observed  at  x  due  to  the 

wave  scattered  by  the  s^*1  inclusion  located  at  x  and  where 
n 

u  (x)  is  the  displacement  due  to  the  incident  wave. 

When  examined  on  the  scale  of  V^,  Figure  1  also  represents 
a  piece  of  a  two-phase  medium  isolated  in  an  infinite  matrix. 
Assuming  that  the  two-phase  medium  is  homogeneous  on  the 
scale  of  Vq  we  can  define  the  properties  of  the  effective 
medium  to  be  the  same  as  those  of  an  homogeneous  medium 
which,  when  confined  to  the  volume  VQ  and  illuminated  by  the 
same  incident  wave  u°  produces  the  same  displacement  field 
at  the  point  x  as  the  field  generated  by  the  N  inclusions. 

We  may  write 

u(x)  ■  u°(x)  +  u* (x  ,x^)  (2) 

where  u*(x  ,Xg)  is  the  scattered  displacement  field  observed 
at  x  due  to  the  volume  Vq  having  effective  properties  and 
located  at  xQ. 

Equating  (1)  and  (2)  we  obtain  the  fundamental  equation 
defining  the  effective  medium 


uMx,)^) 


N 

E  u*(x  ,x  ) 

s-1  ~ 


(3) 


If  we  assume  that  the  two-phase  medium  is  isotropic,  the 


effective  medium  will  also  be  isotropic  and  we  must  take  a 
sphere  for  the  volume  VQ  so  that  the  scattered  waves  do  not 
depend  on  the  orientation  of  VQ  with  respect  to  the  incident 
field.  Thus  from  here  on  we  shall  call  the  volume  VQ  the 
representative  sphere. 

In  order  to  solve  (3)  for  the  effective  properties 
exactly/  a  complete  statistical  description  of  the  distri¬ 
bution  of  the  inclusions  is  required.  The  wave  scattered 
by  each  inclusion  is  a  function  of  the  wave  incident  on 
this  particular  inclusion  and  this  incident  wave  depends  on 
the  location  of  the  other  inclusions  because  of  multiple 
scattering  effects.  The  relative  location  . *  all  inclusions 
must  also  be  known  for  the  calculation  of  the  sum  of  all 
scattered  waves.  In  real  two-phase  media  such  statistical 
information  is  not  available.  Usually  only  the  relative  volume 
fractions  of  the  phases  are  known.  In  order  to  solve  (3)  with 
the  limited  information  we  have,  we  make  two  additional 
assumptions:  (i)  we  assume  that  the  observation  point  x  is 

sufficiently  far  from  the  representative  sphere  so  that  we  may 
take  as  a  first  approximation 

x^  *  /  s  -  1,  N 

and  (ii)  we  assume  that  multiple  scattering  effects  are 
negligible,  which  allows  us  to  take  the  undisturbed  incident 
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field  as  the  field  incident  on  each  inclusion  within  the 
representative  sphere.  The  latter  assumption  restricts  the 
validity  of  our  results  to  small  volume  concentrations  of 
inclusions ,  or  in  other  words  to  the  case  of  non-interacting 
inclusions.  The  possibility  of  taking  interactions  into 
account  by  the  use  of  self-consistent  schemes  is  discussed 
elsewhere  (Kuster,  1972) . 

With  these  approximations,  the  problem  of  finding  the 
effective  properties  of  a  two-phase  medium  reduces  to  the 
estimation  of  the  scattered  displacement  field  due  to  an 
inclusion  isolated  in  an  infinite  matrix,  given  a  mono¬ 
chromatic  incident  wave  of  long  wavelength. 

Even  thif  simplified  problem  is  amenable  to  mathematical 
treatment  only  when  the  inclusion  is  of  a  regular  shape,  such 
as  a  sphere  or  a  spheroid. 


Scattering  by  spherical  inclusions 

Consider  an  homogeneous  and  isotropic  sphere  with 
elastic  constants  X'  and  u'  and  density  p'  embedded  in  an 
homogeneous  and  isotropic  infinite  matrix  with  elastic  con¬ 
stants  X  and  u  and  density  p.  Let  a  plane  P  wave  be  incident 
along  the  x  axis  with  displacement 


i  (px-wt) 

v 


(4) 
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where  A  is  the  amplitude,  p  the  wavenumber  and  w  the  angular 
frequency.  The  presence  of  the  sphere  generates  four  additional 
waves:  the  P  and  S  waves  inside  the  sphere  and  the  P  and  S 
waves  scattered  into  the  matrix.  We  formulate  the  problem 
in  spherical  coordinates  with  the  origin  at  the  center  of  the 
sphere  (Figure  1) .  Because  of  the  symmetries  of  the  problem 
we  need  not  consider  the  azimuthal  dependence.  Following 
Yamakawa  (1962)  we  express  the  radial  and  transverse  dis¬ 
placements  u  and  v  corresponding  to  each  wave  in  infinite 
series  of  spherical  Bessel  functions  and  Legendre  polynomials. 
Thus  the  incident  P  wave  components  are : 


V*4  I 

p  n-0 


0"“S  _L  (2n  *  1)10  ^  jn(pr)Pn(cos0) 


A  7  ,  ,wi>  ^n(pr)  d  „  , 
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The  scattered  P  waves  are 
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The  scattered  S  waves  are 

l  m  n~' 

u2  - - J  2  C  n(n  +  1)  -S_ 

2  n«l  n  r 


hi1'  Ur, 
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The  P  waves  in  the  sphere  are 


u. 


£  0 


p'  n-0 
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The  S  waves  in  the  sphere  are 


u4  *  -  Z  E„  "(n  ♦  1)  -2— 

s'2  n-1  n  r 


P„ (cose) 
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1  T  En  d 
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(9) 


The  e  time  dependence  is  omitted  throughout  for  brevity 
and  p  and  s  ( p *  and  s')  denote  the  wavenumbers  of  P  and  S 
waves  in  matrix  (inclusion)  material.  Pn(cos8)  is  the 
Legendre  polynomial  of  the  nth  order,  jn(z)  is  the  spherical 
Bessel  function  of  the  nth  order  and  hU>  (x)  is  th#  8pherical 

Hankf.l  function  of  the  first  kind  and  the  nth  order.  We  must 
use  hiJ1)  (z)  for  waves  travelling  radially  outward  because 
we  adopted  an  e  iut  time  dependence.  The  coefficients  in  the 
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series  are  determined  from  the  boundary  conditions  on  the 
surface  of  the  sphere  (continuity  of  displacements  and 
normal  stresses  at  r  ■  a) .  The  detailed  calculation  of  the 
scattered  waves  is  given  in  Appendix  A  in  the  case  of 
interest  to  us  vhere  the  wavelengths  of  all  waves  (incident 
as  well  as  scattered  and  transmitted)  are  much  longer  than 
the  radius  of  the  sphere  and  where  the  observation  point  is 
at  a  large  distance  from  the  sphere.  With  these  approximations 
the  scattered  waves  can  be  written  as 
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iA  .  .3  e* (Pr~ut)  Be, 
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The  neglected  terms  are  of  order  (pa)5. 

It  was  shown  by  Yamakawa  (1962)  that  in  the  case  of  a 
spherical  cavity  the  scattered  waves  can  be  obtained  from 
(10)  and  (11)  by  letting  K',  y'  and  p'  vanish; and  in  the 
case  of  fluid-filled  cavity  by  letting  u’  vanish.  However 
the  solution  for  the  case  of  a  fluid  matrix  cannot  be  obtained 
by  simply  letting  y  vanish.  In  this  case,  *-he  scattered  waves 
can  be  written  as  (see  Appendix  A  for  the  detailed  derivation) 
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The  terms  neglected  are  of  order  (pa)  .  This  result 
holds  whether  the  inclusion  is  solid  or  fluid. 

The  difference  between  the  solid  and  fluid  matrix  case 
resides  essentially  in  the  2nd  term  (n  -  1) .  This  term 
represents  a  single  force  source  which  expresses  the  change 
in  inertia  due  to  the  replacement  of  matrix  material  by 
inclusion  material.  When  the  matrix  is  solid  the  change  in 
inertia  arises  only  from  the  density  difference  between 
matrix  and  inclusion  since  there  is  no  relative  motion  between 
the  matrix  and  inclusion.  But  when  the  matrix  is  fluid 
relative  motion  does  occur  (Lamb, 1932  )  and  the  inertia  term 
is  therefore  modified. 

A  spheroidal  inclusion 

We  wish  to  find  the  waves  scattered  by  a  spheroidal 
inclusion  embedded  in  a  solid  matrix.  We  cannot  follow  the 
same  procedure  as  for  a  spherical  inclusion  since  the  vector 
wave  equation  is  not  separable  in  spheroidal  coordinates 
(Morse  and  Feshbach,  1953).  Rather  we  shall  use  an  integral 
expression  derived  by  Mai  and  Knopoff  (1967)  for  the  displacement 
due  to  the  waves  scattered  by  an  inclusion  of  arbitrary  shape 
isolated  in  an  infinite  matrix.  Denoting  the  displacement 
field  of  the  scattered  waves  observed  at  a  point  x  in  the 
matrix  by  Au(x) ,  the  displacement  field  of  the  wave  at  a 
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point  £  inside  the  inclusion  by  v(£) ,  the  kth  component  of 

the  Green's  function  due  to  a  point  force  acting  in  the  ith 

direction  at  a  point  ^  in  the  infinite  matrix  by  Gk^(x,y;) 

and  the  elastic  tensor  in  the  matrix  by  c. .  ,  we  have 

ijpq 


Au^ (x) 


ip'-p)  vi(£)Gki(x,£) 


3v_  3G 


“(cijpq  “  Cijpq# 


-Irr  >  * 


(14) 


The  summation  convention  is  used  and  we  omit  for  brevity 
the  assumed  e*wt  time  dependence.  The  integral  is  taken  over 
the  volume  V  of  the  inclusion.  If  the  displacement  and 
strain  inside  the  inclusion  can  be  estimated  in  terms  of  the 
incident  field  we  can  obtain  the  desired  expression  for  the 
scattered  field. 

The  estimation  of  the  displacement  and  strain  inside 
the  inclusion  is  given  in  Appendix  B  under  the  assumption  that 
the  matrix  is  solid,  that  the  inclusion  is  spheroidal  and  that 
the  wavelengths  are  much  longer  than  the  inclusion  size.  For 
an  observation  point  located  at  a  large  distance  from  the 
inclusion,  the  displacement  field  scattered  by  a  spheroid 
of  arbitrary  orientation  can  be  written  as 
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The  corresponding  expression  for  a  spherical  inclusion 
was  derived  by  Mai  and  Knopoff  (1967)  . 

Auk(x,£)  -  V  j^w2(p'-p)uJ(£)Gki(x,;)  - 
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THE  EFFECTIVE  PROPERTIES  OF  TWO-PHASE  MEDIA 

In  order  to  obtain  the  effective  properties  of  two-phase 
media  we  row  use  equation  (3)  and  the  expressions  derived 
above  for  the  scattered  waves. 

Spherical  inclusions 

Consider  N  spherical  inclusions  embedded  in  an  infinite 
solid  matrix  and  confined  to  the  representative  sphere  of 
radius  R.  Let  a  plane  P  wave  of  amplitude  A  be  incident 
from  infinity.  At  a  distant  observation  point,  with  the 


-63- 


Auk(x,;)  -  v 


2 

w  ( p  *  **p)  u 


i<£>Gki‘*'£> 


0  3Gki 

+  2(m,-m)U.  .  ]e°  -r£±  ( 

H  i  jrs  rs  1 


x,;)J 


(15) 


The  corresponding  expression  for  a  spherical  inclusion 
was  derived  by  Mai  ai  d  Knopoff  (1967)  . 
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THE  EFFECTIVE  PROPERTIES  OF  TWO-PHASE  MEDIA 

In  order  to  obtain  the  effective  properties  of  two-phase 
media  we  now  use  equation  (3)  and  the  cxprrssions  derived 
above  for  the  scattered  waves . 

Spherical  inclusions 

Consider  N  spherical  inclusions  embedded  in  an  infinite 
solid  matrix  and  confined  to  the  representative  sphere  of 
radius  R.  Let  a  plane  P  wave  of  amplitude  A  be  incident 
lrom  infinity.  At  a  distant  observation  point,  with  the 
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assumptions  listed  in  the  previous  section,  the  waves 
scattered  by  the  representative  sphere  are  given  by  an 
expression  of  the  type  of  equation  (10) 

* 
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Assuming  that  multiple  scattering  effects  arc  negligible  so 
that  the  wave  incident  on  each  inclusion  is  the  original 
plane  P  wave  we  can  also  write  the  scattered  waves  at  the 
obsarvation  point  as  the  sum  of  the  wave-,  scattered  by  each 
inclusion.  Thus  we  also  have 
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Because  of  our  definition  of  the  effective  medium  by  equation 
(3) ,  we  obtain  the  effective  properties  by  equating  (17)  and 
(18).  Since  the  equality  must  hold  independently  of  the  angle 
0  the  coefficients  of  the  corresponding  angular  terms  must 
be  equal.  Thus  we  obtain  the  composition  laws  for  the 
effective  clastic  constants  and  density 


K*-K  _  K’-K 
3K*+4y  3Kj  +  4'< 


(19) 


p*-p  -  c(p'-p) 


(20) 


(21) 


where  c,  the  volume  concentration  of  inclusions  is  given  by 

i  ;  3 

c  -  t  1  ai  (22) 

K  j-1  3 

The  corresponding  effective  P  and  S  wave^veloci ties  are 

6*-M? 

The  effective  properties  derived  here  were  obtained  earlier 
by  Ament  (1959) .  The  assumptions  involved  in  their  derivation 
are;  (i)  the  matrix  is  solid,  (ii)  the  inclusions  are 
spherical,  (iii)  the  wavelengths  of  all  waves  arc  much  longer 
than  the  inclusion  radius  and  (iv)  multiple  scattering  effects 
can  be  neglected.  Because  of  the  latter  assumption,  the 
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validity  of  this  model  is  limited  to  two-phase  media  where 
the  concentration  of  inclusions  is  small.  Als-o,  the  effective 
clastic  moduli  in  this  model  reduce  to  those  found  by  Mai 
and  Knopoff  (1967)  if  the  concentration  of  inclusions  is  much 
smaller  than  unity.  The  Mai  and  Knopoff  model  is  probably 
valid  for  lower  concentrations  than  the  model  derived 
above . 

In  the  case  of  a  non-viscous  fluid  matrix  one  must  use 
;12)  as  the  expression  of  the  waves  scattered  both  by  the 
representative  sphere  and  by  each  solid  spherical  inclusion. 
Introducing  (12)  with  the  appropriate  starred  or  primed 
variables  in  (3)  ,  identifying  the  coefficients  of  the 
corresponding  angular  terms  and  using  (22)  we  obtain  the 
following  composition  laws 


K-K*  „  K-K' 
■**“  C 


(24) 


p- P*  D-D ' 

p+lp7  "  c  p73FT  (25) 

The  effective  shear  modulus  vanishes  because  a  two-phase 
medium  cannot  sustain  any  shear  unless  there  is  a  solid 
continuum.  The  effective  P  wave  velocity  is  given  by 

a*  ■  /k*/P*  (26) 

It  is  of  interest  to  note  that  the  effective  P  wave  velocity 
is  independent  of  the  shear  modulus  of  the  inclusions.  Thus 
our  result  is  identical  with  that  obtained  by  Ament  (1953) 

for  suspensions  of  perfectly  rigid  spheres  in  a  non-viscous  fluid 
matrix. 


r 
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It  must  bo  emphasized  that  p*  is  an  effective  inertial 
density  since  (25)  is  a  relation  between  inertia  terms.  In 
the  case  of  a  fluid  matrix  where  relative  motion  between 
the  inclusion  and  mp^rix  can  occur,  the  effective  inertial 
density  is  different  from  the  effective  gravitational  density 
which  is  given  by 


p*  »  p(l-c)  +  p'c 


(27) 


Because  of  this  difference  resulting  from  inertia  effects, 
the  static  and  dynamic  velocity  formulae  vary.  In  the  static 
problem  the  effective  bulk  modulus  is  given  by  the  Reuss 
average  which  is  in  fact  equation  (24),  but  one  would 
have  taken  p*  as  the  density  of  the  medium  and  thus  the 
effective  "static"  P  wave  velocity  would  have  been 


The  "static"*  velocity  calculated  with  equation  (28)  can  be 
significantly  different  from  the  effective  P  wave  velocity 
given  by  equation  (26)  when  the  density  contrast  between  matrix 
and  inclusion  materials  is  large. 

Spheroidal  inclusions 

We  again  use  equation  (3)  to  obtain  the  effective 
properties.  The  field  scattered  by  the  representative  sphere 
centered  at  ^  is  given  by 
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where  P*  and  Q*  are  obtained  from  (B-4)  by  letting  all 
primed  variables  become  starred  variables.  In  order  to 
evaluate  the  sum  of  tho  fields  scattered  by  all  spheroids 
within  the  representative  sphere,  we  assume  that  (i)  multiple 
scattering  effects  are  negligible  so  that  the  field  incident 
on  each  spheroid  is  also  u°,  ( ii )  the  orientation  of  the 
spheroids  is  uniform  over  all  directions  so  that  the  two- 
phase  medium  is  isotropic  on  a  large  scale,  (iii)  the 
distribution  function  of  the  orientation  can  be  represented 
by  a  continuous  function  although  there  is  a  finite  number  of 
sph  sroids  in  the  representative  sphere,  and  (iv)  all 
inclusions  arc  approximately  located  at  £°.  Then  using  (is) 

wo  may  write  the  sum  of  the  fields  scattered  by  all  inclusions 

as 
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where  is  the  volume  of  the  nth  inclusion,  N  is  the  number 
of  inclusions  in  the  representative  sphere  and 
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Equating  (29)  and  (30)  and  setting 
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Rcquiring  that  both  terms  are  identically  zero,  we  have  the 
density  composition  law  from  the  first  term 


p*  ■  p(l-c)  +  p'c 


(34) 


and  the  composition  law  for  the  elastic  constants  from  the 
second  term.  If  the  incident  field  is  purely  dilatational , 
that  is 
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the  second  term  becomes 
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From  the  symmetries  of  the  tensor  (see  Appendix  3 

and  the  integration  in  (31)  we  find  that 
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Combining  (36)  and  (37)  we  obtain  the  composition  law  for 
the  bulk  modulus 
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The  composition  law  for  the  shear  modulus  is  obtained  in  a 
very  similar  way  by  taking  the  incident  field  as 
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Thus  wc  have 
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The  scalars  T^^  and  T^^  are  functions  of  the  shape  of  the 
spheroid.  They  are  given  in  tho  Appendix. 

The  essential  result  here  is  that  the  effective  elastic 
moduli  depend  not  only  on  the  concentration  but  also  on  the 
shape  of  the  inclusions  (i.e.  aspect  ratio).  This  conclusion 
is  in  agreement  with  experimental  data  on  porous  rocks 
(Nur  and  Simmons,  1969)  and  with  other  theoretical  formulations 
(Eshelby,  1957;  Wu,  1966;  Walsh,  1969).  It  is  important 
to  note  that  our  assumption  of  non-interaction  among  the 
spheroids  is  violated  when  the  ratio  c/a  is  larger  than  1 
(Solomon,  1971)  since  the  inclusions  are  then  overlapping, 
at  least  partially. 


In  the  above  formulation,  the  results  were  for 
spheroids  having  all  the  same  aspect  ratio.  These  results 
can  be  easily  extended  to  cover  the  case  of  a  discrete 
spectrum  of  aspect  ratios.  When  the  concentration  for  each 
aspect  ratio,  c(am)  is  known,  the  effective  bulk  modulus 
is  given  by 


K*-K 


K  1  -K 
3K+Tu 


M 

E 

m«l 


c(a  ) 
m 


1 

3 


Tii j j 


(41) 


-72- 


0£  course,  the  non-interaction  assumption  must  still  be 
valid  and  it  can  be  expressed  as 
M  c(a  ) 

E  - E-  <  l  ( 

m»l  am 

We  can  compare  our  results  (38)  and  (40)  with  those 
obtained  by  Walsh  (1969)  for  a  two-phase  medium  with  non 
interacting  spheroidal  inclusions. 

K*-K  "  I  Tiijj(K'"K) 

x  < 

y*-y  -  |  (y'-M) <Ti ji j  "  J  Tiijj) 

To  compare  these  with  our  results ,  we  rewrite  (38)  and 
(40)  in  a  slightly  different  form, 


K*-K 


C 

I 


(K,-K)Tiijj 


3K*+4u 

3K+ly 


(44) 


I  ^•-y)CTijlj 


iiil.  6u* (Kt2u)tu(9Kt8u) 
~y^}  5u(3Kt4u) 


It  is  clear  that  our  results  and  those  of  W^lsh  differ 
somewhat.  They  are  similar  if  the  effective  medium  and  the 
matrix  arc  not  too  different.  This  may  be  the  case  if  the 
concentration  of  inclusions  and/or  the  clastic  moduli 
contrast  between  inclusion  and  matrix  materials  is  small, 
if  the  effective  moduli  derived  for  spheroidal  inclusions 
are  specialized  to  the  limiting  case  of  spheres  where 


Now 
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1  t  .  3K4- 4m 
3  Aiijj  3K’+4y 


2Su  (3K+4p) 

6m  1  (K+2ii)+vi(9K+8u) 


(45) 


the  results  of  Walsh  reduce  to  those  of  Mai  and  Knopoff 
whereas  our  results  (i.e.(44)  become  (19)  and  (21).  Thus 
at  least  in  the  case  of  spherical  inclusions,  our  results 
are  probably  valid  over  a  wider  range  of  concentration  than 
those  of  Walsh  (1969)  . 


NUMERICAL  RESULTS  AND  DISCUSSION 

To  demonstrate  the  effects  of  the  inclusions  on  th*: 
effective  moduli  and  velocities  of  a  two-phase  medium  w  i 
made  a  series  of  numerical  calculations.  For  these 
calculations  we  adopted  a  solid  matrix  and  inclusions  of 
different  shapes  and  properties  to  represent  both  dry  and 
water-saturated  rocks.  The  matrix  parameters  were  chosen 
to  represent  those  of  the  matrix  of  an  average  sandstone 
or  quartz  rich  crystalline  rock:  K  -  0.44  Mb,  u  ■  0.37  Mb, 
p  ■  2.70  g/cc.  For  the  "saturated"  case  we  used  water 
inclusions  with  the  moduli:  k'  -  0.022  Mb,  u  “0, 
p'  ■  1.00  g/cc.  In  the  case  of  the  "dry"  state,  it  was 
assumed  that  the  pores  were  filled  with  air  at  atmospheric 
pressure  (bulk  modulus  k'  -1.5  bars) . 
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With  these  input  parameters,  two  separate  typer  of 
calculations  were  carried  out.  The  first  set  of  calculations 
was  to  show  the  effect  of  each  pore  type.  For  this  we 
specified  the  pore  shape  and  calculated  the  moduli  and 
velocities  of  the  composite  medium  as  a  function  of 
concentration  of  pores  for  both  the  "saturated"  and  "dry" 
states.  The  results  are  shown  in  Figures  2a, b. 

The  four  aspect  ratios  chosen  are  a  ■  1.0  (spherical 
—1  -2  -4 

pores)  and  a  >•  10  ,10  ,10  ,  representing  oblate 

spheroidal  pores  of  different  shapes  going  all  the  way  into 
fine  cracks.  For  ellipsoidal  pores  the  concentrations  were 
varied  from  c  ■  0  to  c  ■  a,  the  optimum  limit  specified 
by  non-interaction  assumption  (equation  (4.:)  1  .  The 
densities  were  computed  using  equation  (27) . 

The  four  examples  shown  in  Figures  2a, b  demonstrate 
the  effects  of  inclusion  concentrations  and  shapes  on  the 
effective  moduli  and  the  velocities.  In  all  cases  the 
moduli  (K* ,  y*)  and  the  velocities  decrease  with  increasing 
concentration  of  inclusions.  For  a  given  concentration  the 
flatter  inclusions  have  relatively  greater  effect  than  the 
rounder  inclusions.  Even  a  very  low  concentration  of  thin 
inclusions  (i.e.  0.01  percent)  could  decrease  the  velocities 
in  the  composite  medium  by  as  much  as  ten  percent  or  more. 

The  comparison  of  the  moduli  for  the  water- saturated  or 


"dry"  states  reveals  some  important  effects.  At  a  given 
concentration  of  pores,  when  we  go  from  water  saturated  to 
the  air  saturated  ("dry")  state  the  relative  change  in  bulk 
modulus  is  greater  than  the  corresponding  change  in  shear 
modulus.  This  is  true  regardless  of  the  pore  shape.  For 
flat  pores,  the  changes  in  both  moduli  are  pronounced  while 
for  spherical  pores  only  the  change  in  the  bulk  modulus  is 
apparent. 

The  effects  of  water  or  gas  saturation  on  compressional 
and  shear  velocities  strongly  depend  on  pore  shapes  as 
illustrated  in  Figures  2a, b.  For  spherical  pores  the 
velocities  increase  when  one  goes  from  the  water  saturated 
to  rary!  state.  This  is  because  the  effect  of  density  change 
in  the  composite  medium  is  greater  than  the  changes  of  the 
moduli.  Thus  velocities  increase  in  dry  state  while  bulk 
and  shear  moduli  are  decreasing,  because  of  the  greater 
decrease  in  effective  densities  (equation  23)  .  For  inter¬ 
mediate  shape  pores  such  as  those  with  aspect  ratio  of  about 
a  ■  0.1,  the  velocities  in  the  "dry"  and  water  saturated  states 
are  about  the  same.  For  very  thin  pores  and  cracks,  both  P 
and  S-velocitios  are  lower  in  the  "dry"  state  than  in  the 
water-saturated’  case,  although  the  decrease  is  much  more 
pronounced  in  the  case  of  P-waves.  In  these  cases  tht 
density  changes  are  negligibly  small;  velocities  in  two  stages 
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are  directly  controlled  by  the  effective  moduli. 

In  a  typical  rock  the  pores  are  likely  to  represent  a 
spectrum  of  shapes.  Equidimensional  pores  in  sedimentary 
rocks  can  be  approximated  by  spheres  while  grain  boundary 
spaces  can  be  approximated  by  flat  cracks  (low  aspect  ratio 
spheroids)  as  shown  by  electron  microscope  studies  (Timur 
et  al.,  1972).  To  study  the  effects  of  saturation  in 
rocks,  a  second  set  of  calculations  was  carried  out  for 
models  of  a  sedimentary  rock  (sandstone)  and  a  crystalline 
rock  with  representative  porosities:  14.2  percent  for 
sedimentary  and  0.4  percent  for  crystalline.  The  matrix 
model 
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and  the  properties  of  the  fluids  (water  and  air)  were  the 
same  as  before .  A  spectrum  of  pore  shapes  was  u&od  in  each 
case.  For  the  sedimentary  model  the  majority  of  pores  (12 
percent  concentration)  were  assumed  to  be  spherical  pores; 
the  remainder  had  smaller  aspect  ratios.  These  are  listed 
in  Table  1.  For  the  crystalline  model  the  majority  of  the 
pores  were  taken  to  be  thin  cracks  as  indicated  by  the  thin 
section  photographs  (Timur  et  al.,  1972;  Brace  et  al., 

1972)  . 

The  calculated  moduli  and  velocities  are  listed  in 
Table  1.  In  all  cases  calculations  were  carried  out  for 
the  standard  1  atmosphere  pressure.  The  differences  in 
moduli  and  velocities  between  the  water  saturated  and  "dry'' 
states  arc  very  significant.  For  the  sedimentary  rock 
model  the  compressional  velocity  increases  about  22  percent, 
relative  to  the  "dry"  state  wnen  saturated  with  water.  The 
shear  velocity  changes  about  3  percent.  For  the  crystalline 
rock  model,  althjugh  the  total  porosity  is  only  0.4  percent,  the 
increase  of  ccmpressional  velocity  upon  saturation  is  about 
15  percent  while  the  change  in  shear  velocity  is  about  4 
percent.  Velocity  changes  of  this  nature  have  been  observed 
in  the  laboratory  for  dry  and  saturated  granites  (Nur  and 
Simmons,  1969). 

Those  two  sets  of  examples  demonstrate  the  relative 
importance  of  the  pore  shapes  and  the  compressibility  and 
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and  density  of  saturating  fluids  on  determining  the  velocities 
in  rocks.  If  pore  shapes  can  be  specified  on  the  basts  of 
laboratory  measurements,  then  the  nature  of  saturating 
fluids  could  be  determined  from  P  and  S-wave  velocities. 

The  attenuation  of  clastic  waves  in  two-phase  media 
can  be  computed  using  the  formulations  given  in  this 
paper  by  assuming  the  moduli  (K,  u)  are  complex.  This  is 
discussed  in  detail  in  the  second  paper  (Kuster  and  Toksdz , 
1973b)  along  with  experimental  results. 


APPENDIX  A 


Consider  an  homogeneous  and  isotropic  sphere  with 
elastic  constants  X'  and  p  '  and  density  p'  embedded  in  an 
homogeneous  and  isotropic  infinite  matrix  with  elastic 
constants  X  and  p  and  density  p.  Let  a  plane  P  wave  bo 
incident  along  the  x  axis  with  displacement 
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whore  A  is  the  amplitude,  p  the  wavenunber  and  w  the  angular 
frequency.  The  presence  of  the  sphere  generates  four  additional 
waves:  the  P  and  S  waves  inside  the  sphere  and  the  P  and  S 

waves  scattered  into  the  matrix.  We  formulate  the  problem 
in  spherical  coordinates  with  the  origin  at  the  center  of  the 
sphere  (Figure  1) .  Decausc  of  the  symmetries  of  the  problem 
we  need  not  consider  the  azimuthal  dependence.  Following 
Yamakawa  (1962)  we  express  the  radial  and  transverse  dis¬ 
placements  u  and  v  corresponding  to  each  wave  in  infinite 
series  of  spherical  Desscl  functions  and  Legendre  polynomials. 
Thus  the  incident  P  wave  components  arc: 
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The  scattered  P  waves  are 
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The  scattered  S  waves  arc 
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The  P  waves  in  the  sphere  are 
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The  S  waves  in  the  sphere  are 
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(A- 6) 


The  e”iwt  time  dependence  is  omitted  throughout  for  brevity 

and  p  and  s  ( p '  and  s')  denote  the  wavenumbers  of  P  and  S 

waves  in  matrix  (inclusion)  material.  Pn(cos0)  is  the 

Legendre  polynomial  of  the  n^  order,  3n(z)  spherical 

Bessel  function  of  the  nth  order  and  h^1*  (2)  is  the  spherical 

Hankel  function  of  the  first  kind  and  the  nth  order.  We  must 

use  h^(z)  for  waves  travelling  radially  outward  because 
n 

we  adopted  an  o”^ut  time  dependence.  The  coefficients  in  the 


series  are  determined  from  the  boundary  conditions  on  the 

surface  of  the  sphere  (r  ■  a) .  The  boundary  conditions 

are  continuity  of  the  displacements  and  of  the  normal  stresses 
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where 
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Putting  (A-2)  to  (A-6)  in  (A-7)  we  obtain  the  following  system  of 
four  equations  with  four  unknowns  for  each  nil 

BnlB„  +  YnlCn  +  6nlDn  +  C„1E„  ’  , 

Bn2B„  +  Yn2Cn  *  4n2Dn  +  cn2E„  ’  i’Wlla^A  , 

Bn3Bn  *  Yn3Cn  +  in3t>„  +  cn3En  ’  i" <*■•+!> °n3A  - 


>  (A-8) 
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For  n  -  0,  equation  (a-8)  can  be  written  as 
601B0  +  601D0  “  a01A 
603B0  +  603D0  “  a03A 

where 

a0l  "  ”  l  ( X+2u)  (O  - 
«C:  *  l(»+2g)h^1,(C)  -  ^  h‘l,(t)l 
<01  "  -t(X>+2u')30(C‘)  -  ^U1)) 

and 

“03  '  -  l 
«03  *  X  h{1’(«' 

<03  '  -  b  hiV) 


(A-13) 


(A-14) 


(A-15) 


Cramer's  rule  can  be  used  to  solve  systems  (A-8)  and  (A-13) 
for  the  coefficients  in  the  series  expansions  of  the  scattered 
waves.  When  the  wavelengths  of  all  wavcc  are  much  longer  than 
the  radius  of  the  sphere  (C,  n,  n '  much  smaller  than 
1)  wo  may  use  the  expansions  of  the  spherical  Bessel  and 
Hankel  functions  for  small  arguments.  They  are  for  z  <<  1 
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Using  the  expansions  in  (A-9)  to  (A-12) , 

(A-14 )  and  (A-15) 

keeping  only  the  dominant  term  in  solving 

systems  (A-8)  and  (, 

we  obtain 
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For  an  observation  point  at  a  large  distance  from  the  sphere 
(pr  >>  1,  sr  >>  1)  we  may  use  the  asymptotic  expansion  of  the 
Hankcl  function 
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and  we  may  thus  writs  the  scattered  P  waves  as 
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and  the  scattered  S  waves  is 
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Since  pr  >>  1  and  sr  >>  1,  v^  car.  be  neglected  w*  th  respect 
to  u^,  and  u2  can  be  neglected  with  respect  to  v2-  Further¬ 
more  in  the  long  wavelength  case  all  coefficients  for  n  1  3 
c*e  of  higher  order  than  the  leading  coefficients.  Thus, 
keoping  only  the  dominant  terms,  we  can  write  the  scattered 


waves  as 
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When  the  matrix  is  a  non-viscous  fluid,  there  are  no 
scattered  S  waves  and  the  continuity  of  the  transverse  dis¬ 
placement  at  the  boundary  is  not  required.  Thus  the  coefficients 
in  the  expansion  of  the  scattered  P  waves  arc  the  solutions 
of  the  following  systems. 


4 
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for  n  ■  0 
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for  nil 
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The  coefficients  6^  and  are  the  same  a3  in  (A-9)  to 

(A-ll)  whereas  6®  and  are  obtained  from  B  and  a  in 

nm  nm  nm  nm 

(A-9)  to  (A-ll)  by  letting  u  vanish.  The  calculation  is 
now  similar  to  that  of  the  preceding  section.  In  the  long 
wavelength  approximation  ($,  C»  and  n'  much  smaller  than  1) 
and  for  an  observation  point  at  a  large  distance  from  the 
sphere  (pr  >>  1) ,  the  scattered  waves  can  be  written  as 
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The  terms  neglected  are  of  order  (pa)5. 
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APPENDIX  B 

To  find  the  waves  scattered  by  a  spheroidal  inclusion 
embedded  in  a  solid  matrix,  we  shall  use  an  integral 
expression  derived  by  Ma'  and  Knopoff  (1967)  for  the 
displacement  due  to  the  waves  scattered  by  an  inclusion  of 
arbitrary  shape  isolated  in  an  infinite  matrix.  Denoting 
the  displacement  field  of  the  scattered  waves  observed  at 

a  point  x  in  the  matrix  by  Au(x)  ,  the  displacement  field 
of  the  wave  at  a 
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point  £  inside  the  inclusion  by  v(£)  ,  the  kth  component  of 
the  Green's  function  due  to  a  point  force  acting  in  the  ith 
direction  at  a  point  £  in  the  infinite  matrix  by  G^(x,y) 
and  the  elastic  tensor  in  the  matrix  by  c^  ,  we  have 


(B-l) 


The  summation  convention  is  used  and  we  omit  for  brevity 
the  assumed  e^ut  time  dependence.  The  integral  is  taken  over 
the  volume  V  of  the  inclusion.  If  the  displacement  and 
strain  inside  the  inclusion  can  be  estimated  in  terms  of  the 
incident  field  we  can  obtain  the  desired  expression  for  the 
scattered  field. 

When  the  matrix  is  solid  and  when  the  wavelengths  are 
much  longer  than  the  inclusion  size,  the  lowest  order 
approximation  (known  as  the  Born  Approximation)  to  the  dis¬ 
placement  inside  is  the  displacement  one  would  observe  if 
the  inclusion  were  absent  (Mai  and  Knopoff ,  1967)  .  In 
other  words,  we  may  write 

v(£)  =  u°  (O  (B-2) 


where  £  is  the  center  of  the  inclusion.  No  assumption  about 
the  shape  of  the  inclusion  is  required  so  that  we  can  use 
(B-2)  for  a  spheroidal  inclusion. 
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Mal  and  Knopoff  (1967)  showed  that  for  arbitrary 
contrast  between  the  elastic  properties,  the  lowest  order 
approximation  to  the  strain  inside  a  spherical  inclusion  is 
terms  of  the  incident  strain  is  given  by 

eklU)  “  j(p"Q)cii 


where  6.  is  a  Kronecker  delta  and 

I*  M 


P 


3K+4p 


(B-4) 


5u (3K+4y) _ 

Q  "  6 W  '  (K+2p)+y  (9K+6wy 


The  important  feature  of  this  result  is  that  it  is  the 
same  as  that  obtained  by  Eshelby  (1957)  for  the  stati--  strain 
inside  a  sphere  when  a  uniform  strain  is  applied  at  infinity. 

Thus  for  a  spherical  inclusion  and  for  waves  of  long  wavelengths, 
the  lowest  order  approximation  is  given  by  the  solution  of 
the  corresponding  static  problem.  We  assume  this  identity 
also  holds  for  a  spheroidal  inclusion.  The  expression  derived 
by  Eshelby  (1957)  for  the  strain  inside  a  spheroid  of 
arbitrary  orientation  with  respect  to  the  fixed  coordinate 
system  of  the  matrix  is 

0 

®ij  "  Uijki 


(B-5) 
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where 


Uijkl  “  loi  lBj  lyk  ToBy« 

The  4.  are  direction  cosines  and  T  -  .  is  a  fourth 
mn  aPY® 

order  tensor  which  is  described  later  in  this  Appendix. 


Introducin9  (B-2)  and  (B-5)  in  (B-l)  and  using 


Cijpq  ‘  X6ij  ‘pq  +  ,‘(<ip4jq  +  SiqV 


(B-6) 


we  obtain 


Auk(x,i)  -  u2(p'-p)uj<£)  f  Gki(x,&dl 


[(X'-X)U  ra«i^  ♦  2<P,-M)0ijt,l c“J 

V  3 

The  integrals  of  the  Green’s  function  over  the  volume 
of  the  spheroid  can  be  evaluated  easily  when  the  point  of 
observation  is  at  a  large  distance  from  the  spheroid.  With 
this  condition  we  have 


(B-7) 


J  Gki*-'-)d-  *  VG).i(-'^ 


(B-8) 


Finally  by  using  (B-8)  and  the  symmetry  of  the  Green's  function 
we  can  write  the  displacement  field  scattered  by  a  spheroid  of 
arbitrary  orientation  as 
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Auk(x,£)  .  y 


u 


i<£)Gki(x,0 


-t(X'-X)U  f 

PPr»  ij 


(B-9) 


3G. 


*  2<n'-w>u<,  )e° 

i jrs  rs  3^ 


Ta6v6  *s  a  fourth  order  tensor  whose  symmetries  are  for  an 
oblate  spheroid  with  aspect  ratio  o 


Tllll  “  T2222  TH33  *  T2233 

T1122  T2211  T3322  “  T3311 

T1212  "  T1221  “  T2121  “  T2112  (B-10) 

T1313  "  T1331  "  T3113  "  T3131 
T2323  “  T2332  "  T3223  "  T3232 

We  also  have  the  relation 

Tllll  T1122  “  2T1212  *  0  (B-ll) 

The  scalars  and  which  are  used  in  the  text  are 

given  by 
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iijj  F, 


ijij  "  T  Tiijj 


2.1.  F4F5*F6F7~F8F9 
F3  F4  F2F4 


(B-12) 


where 


F^  ■  1  +  A|j(g+$)  -  R(jg  +  §4  -  yj] 

F2  -  1  +  A|l  +  §(g+4>)  -  f(3g+5$)]  +  B{3-4R 
+  |(A+3B) (3-4R)[g  ♦  ♦  -  R(g-$+2$2)] 


r  2 

f3  -  1  +  j[r(2-$)  +  g (R 


->] 


F4  -  1  +  jL  +  g  -  R(g-4>)j 
F5  -  A|R(g+«-j)  -  g]  +  B$(3-4R) 

Fg  «  1  +  A[l  +  g  -  R(g+$)j  +  B(l-$)  (3-- 


'7  '  2  +  t[! 


9$  +  3g  -  R(5$+3g)  +  B$(3-4R) 


F0  -  A^l  -  2R  +  |(R-1)  +  |{5R-3)J  +  B(l-$)(3-< 
F9  -  A  Jg  (R-l)  -  R<|  +  B$(3-4R) 
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TABLE  1 

Theoretical  elastic  moduli  and  seismic  velocities  for  examples 
of  a  sedimentary  and  a  crystalline  rock  model  in  "dry"  and 
"water-saturated"  state.  (In  "dry"  state  the  rock  is 
saturated  with  airO 


PORE  STRUCTURE  USED  FOR  CALCULATIONS 


Sedimentary  model 

Crystalline  model 

Concentration  (%)  Aspect  ratio  Concentration 

12  1  (sphere)  0.01 

2  lO"1  0.15 

0.2  10-2  0.20 

0.2  x  10-1  10"3  0.05 

0.1  x  10-2  10"4  0.10  x  10“3 

(%)  Aspect  ratio 

1 

10*1 

10*2 

10-3 

10*5 

CALCULATED  MODULI  AND 

VELOCITIES 

K  (Mb) 

V  (Mb) 

p (g/cc)  V 

(km/sec) 

P 

V  (km/sec) 
s 

Matrix  Properties 

0.440 

0.370 

2.70 

5.88 

3.70 

Sedimentary  Model 

Water- saturated 

0.276 

0.186 

2.44 

4.62 

2.75 

"Dry"  (air-saturated) 

0.107 

0.163 

2.30 

3.76 

2.67 

Crystalline  Model 

Water  saturated 

0.420 

0.303 

2.693 

5.53 

3.35 

"Dry"  (air-saturated) 

0.260 

0.280 

2.690 

4.85 

3.23 
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FIGURE  CAPTIONS 

Fig.  1  Schematic  diagram  illustrating  the  scattering  of 
a  plane  wave  (u°)  by  a  representative  sphere 
(dashed  circle)  of  the  effective  medium.  Individual 
inclusions  are  outlined  by  solid  lines.  The  scattered 
fields  are  evaluated  at  point  x. 

Fig.  2a, b  Elastic  moduli  (K,u)  and  the  compressional  <P)  and 

shear  (S)  velocities  in  a  composite  medium  with  water 
and  air  filled  pores  as  a  function  of  volume  con¬ 
centration  of  inclusions.  The  four  sets  of 
diagrams  demonstrate  the  effects  of  the  inclusion 
shapes  specified  by  the  aspect  ratio,  a  =  1.0 
corresponds  to  spherical  inclusions;  a  =  10-^ 
represents  very  flat  oblate  spheroidal  inclusions 
modeling  fine  cracks  or  grain  boundaries.  The 
solid  curves  are  the  water-saturated  and  the 
dashed  curves  are  the  air-.cdturated  ("dry")., 
cases.  Note  the  multiplier  in  the  concentration 
scale.  The  matrix  moduli  and  velocities  can  be  read 
at  c  «=  0.  Matrix  density  p  »  2.70  g/cc.  The 
inclusion  moduli  are:  water  -  K'  =  22  kb, 

P'  *  0,  and  p'  =  1.00  g/cc;  and  air  -  K'  =  1.5  bars, 

P '  *  0  and  p '  = 


0.00. 


Aspect  Ratio  *  10"' 


ro  c\J  “  Q  O  Q  O 

O  d  C>  O  <£>  rO 

(q^)  snjnpoiAj  (Das/ixo] )  Xpoo|9A 
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4 .  EARTH  STRUCTURE 

4*1  Structure  and  Tectonics  of  the  Arctic  Region  by 

Edward  Dewey  Chapman  (abstract) 

A  Rayleigh  wave  study  of  the  Arctic  Ocean  has  given 
information  on  the  effect  of  propagation  of  these  waves 
across  the  Arctic.  These  waves  propagate  with  three  dis¬ 
tinct  velocity  functions  across  the  three  regions  of  the 
continental  shelf,  mid-ocean  ridge,  and  deep  ocean  basin. 
The  basin  is  interpreted  to  have  a  sediment  layer  of  2 
kilometers  and  a  crustal  basement  of  15  kilometers  over¬ 
laying  an  oceanic  mantle.  The  Nansen  Cordillera  is  inter¬ 
preted  to  be  a  typical  spreading  center  on  the  basis  of 
these  surface  wave  observations. 

The  border  between  the  North  America  Plate  and  the 
Eurasia  Plate  in  plate  tectonic  theory  is  delineated  by 
examination  of  historical  seismicity  and  regional  geology. 
The  border  goes  from  the  Lena  river  delta,  through  the 
Moma  basin  in  Siberia  to  the  city  of  Magadan.  From  there 
it  continues  down  the  center  of  Sakhalin  and  Hokkaido 
islands  to  a  triple  junction  in  the  Japan  arc.  A  pole  of 
rotation  between  these  plates  is  located  at  approximately 
67°N,  155°E.  This  pole  location  explains  the  current 
seismicity  of  Siberia  and  Sakhalin  island. 
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4*2  P££specting  for  Dead  Slabs  by  Sean  C.  Solomon  and 

Rhett  G.  Butler  (abstract) 

At  subduction  zones  that  have  only  recently  ceased  to 
be  active,  the  lithospheric  slab  may  retain  a  seismic  velo¬ 
city  greater  than  that  of  the  surrounding  mantle  even  after 
the  slab  becomes  seismically  dead.  To  seek  the  subduction 
zone  thought  to  have  been  recently  active  along  the  western 
margin  of  North  America,  we  examined  the  variation  with  pro¬ 
pagation  direction  of  P-wave  travel  time  residuals  from 
sources  at  various  distances  and  azimuths  to  seismograph 
stations  in  Washington  and  California.  The  uncertainty  in 
source  location  and  origin  time  was  removed  by  referring  the 
fravel-time  delay  to  a  nearby  station  overlying  presumably 
more  uniform  mantle.  An  eastward-dipping  band  of  anomalously 
early  arrivals  at  several  stations  on  the  western  flank  of 
the  Sierra  Nevada  and  California  Cascades  may  imply  that  a 
dead  slab  is  present  beneath  northern  California,  though  a 
definitive  conclusion  is  premature  r.t  present  because  of  a 
paucity  of  seismic  sources  in  eastern  North  America.  The 
position  of  the  dead  slab  speculatively  suggested  by  the 
travel-time  data  is  roughly  consistent  with  that  predicted 
by  others  on  the  basis  of  heat  flow  and  geochemistry  in 
the  Sierra  Nevada,  and  the  southward  decrease  in  the  magni¬ 
tude  of  the  travel-time  advance  associated  with  such  a 
slab  is  in  agreement  with  the  history  of  subduction  of  the 
Farallon  plate  as  reconstructed  from  ocean  floor  magnetic 
anomalies  and  continental  tectonic  activity. 
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4 . 3  Seismic  Constraints  on  Ocean-Ridge  Mantle  Structure; 

Anomalous  Fault-Plane  Solutions  From  First  Motions 

by  Sean  C.  Solomon  and  Bruce  R.  Julian 
Surma ry 

Fault-plane  solutions  from  P-wave  first  motion  studies 
of  earthquakes  on  the  crests  of  actively  spreading  mid-ocean 
ridges  often  appear  to  have  non-orthogonal  nodal  planes. 

This  anomaly,  together  with  independent  data  from  surface 
waves  and  shear  waves  from  such  events,  is  most  simply  ex¬ 
plained  as  an  effect  of  propagation  of  the  P  waves  through 
the  laterally  heterogeneous  mantle  beneath  the  ridge.  P-wave 
velocity  models  for  the  ocean-ridge  mantle  are  derived  from 
temperature  models,  the  phase  diagram  for  wet  peridotite,  and 
the  expected  dependence  of  velocity  on  temperature  and  phase. 
Models  that  successfully  account  for  the  nodal  plane  non¬ 
orthogonality  are  characterized  by  pronounced  lateral  tempera 
ture  gradients  and  extensive  partial  melting  at  temperatures 
in  excess  of  the  anhydrous  peridotite  solidus,  in  agreement 
with  models  for  the  petrogenesis  of  tholeiitic  basalts  and 
with  evidence  for  a  region  of  very  low  Q  beneath  ridge  crests 
The  depth  to  the  top  of  the  partially  molten  region  is  not 
well  constrained,  but  the  half-width  of  the  zone  of  extensive 
melting  in  the  north  Atlantic  must  be  several  tens  of  kilo¬ 
meters.  A  prediction  of  these  calculations  is  a  variation  of 
the  apparent  non-orthogonality  of  nodal  planes  with  spreading 
rate  for  ridge-crest  earthquakes;  the  limited  data  are 
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consistcnt  with  this  variation.  The  most  promising  tests 
of  proposed  velocity  models  include  long  seismic  refraction 
profiles  and  studies  of  S-wavo  polarization.  Bf/dy-wave 
travel-time  residuals  observed  teleseismically  i.-om  ocean- 
ridge  earthquakes  are  predicted  to  be  nearly  independent 
of  distance  and  azimuth  and  thus  measurable  only  for  events 
with  independently  determined  origin  times. 
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1.  Introduction 

There  is  a  great  deal  of  current  interest  in  the 
structure  of  the  upper  mantle  beneath  mid-ocean  ridges. 

Of  special  importance  are  inferences  that  can  be  made 
about  mantle  composition  and  temperature  from  observations 
of  phase  boundaries  and  lateral  gradients  of  physical 
properties.  The  temperature  field  is  in  turn  coupled  to 
the  flow  pattern  in  the  mantle,  in  particular  to  the  radial 
extent  of  the  zone  of  upwelling  and  to  the  vertical  velocity 
of  flow  beneath  the  ridge  axis.  In  this  paper  we  show  that 
a  reasonable  interpretation  of  several  peculiar  features  of 
ridge-crest  earthquakes  allows  some  new  constraints  to  be 
placed  on  the  seismic  '-elocity  and  temperature  fields  in  the 
upper  mantle  beneath  spreading  centers. 

Normal-faulting  earthquakes  on  the  crests  of  mid-ocean 
ridges  show  an  unusual  characteristic  which  has  been  noted 
only  in  passing  in  the  literature:  Well-constrained  fault- 
plane  solutions  derived  from  P-wave  first  motions  indicate 
an  apparent  non-orthogonality  of  the  nodal  planes  when 
standard  projections  are  used  to  map  observations  back  to 
the  focal  hemisphere.  We  suggest  in  this  paper  that  this 
anomaly  may  be  ascribed  to  body-wave  path  effects  near  the 
source,  specifically  to  changes  in  propagation  direction 
and  ray  parameter  of  P  waves  caused  by  travelling  through 
a  laterally  heterogeneous  velocity  field  in  the  upper  mantle 
beneath  the  ridge  axis.  By  means 
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of  three-dimensional  ray  tracing  for  two-dimensional  velocity 
models  of  the  ocean-ridge  mantle,  we  calculate  the  corrections 
to  standard  techniques  necessary  to  map  accurately  body-wave 
observations  onto  the  focal  sphere.  Applying  these  correc¬ 
tions  to  data  from  ridge  crest  earthquakes,  we  show  that  the 
P-wave  first  motions  can  then  be  fit  by  two  orthogonal  nodal 
planes, as  expected  for  a  double-couple  source. 

There  are  several  independent  sorts  of  information  that 
can  be  brought  to  bear  on  the  modeling  of  seismic  velocity 
for  the  ocean-ridge  mantle.  Among  these  are  the  temperature 
fields  derived  from  models  for  mantle  flow  beneath  a  spreading 
center,  velocity-temperature  relationships  for  materials 
likely  to  constitute  the  upper  mantle,  and  recent  evidence 
for  a  shallow  zone  of  very  low  Q  localized  beneath  the  ridge 
axis.  Though  a  unique  two-dimensional  velocity  model  cannot  yet 
be  selected,  a  combined  application  of  these  constraints  indi¬ 
cates  that  pronounced  lateral  velocity  gradients,  due  to 
temperature  gradients  and  to  regions  of  partial  melting,  are 
required. 

2.  Focal  Mechanisms  of  Ridge-Crest  Earthquakes 

The  source  mechanisms  of  earthquakes  on  the  axes  of 
actively  spreading  mid-ocean  ridges  have  been  studied  exten¬ 
sively  by  Sykes  (1967,  1970a)  and  his  co-workers,  using 
primarily  the  first  motions  of  P  waves.  Most  such  mechanisms 
can  be  clearly  classified  as  normal  faulting,  with  the 
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inf  errc  d  axis  of  least  compressive  pre-stress  nearly 

horizontal  and  roughly  perpendicular  to  the  strike  of  the 
ridge. 

All  of  these  earthquakes  for  which  both  of  the  nodal 
planes  are  well  known  have  the  anomalous  characteristic 
that  the  nodal  planes  fit  to  firsu-motion  data  projected  with 
a  standard  earth  model  onto  the  lower  focal  hemisphere  are 
not  orthogonal.  Rather  the  downward  dilatational  quadrant 
subtends  a  solid  angle  of  less  than  tt  steradians;  the  angle 
between  the  best  fitting  nodal  planes  is  commonly  60®  to  70®. 
At  least  six  such  mechanisms  (see  Fig.  1  and  Table  1)  have 
been  noted  (Sykes  1967,  1970a;  Thatcher  &  Brune  1971;  Conant 
1972)  ,  although  for  tw.  of  these  solutions  one  of  the  nodal 
planes  was  determined  principally  by  its  required  proximity 
to  stations  recording  small  first  arrivals,  surely  a  risky 
procedure  when  dealing  with  waves  which  must  propagate  beneath 
an  actively  spreading  ridge. 

The  non-orthogonality  of  the  nodal  planes  is  most  likely 
an  artifact  of  the  projection  used  to  map  the  earth's  surface 
back  to  the  focal  sphere.  Specifically  the  standard  tables 
of  ray  parameter,  or  angle  of  incidence,  versus  epicentral 
distance  are  inadequate  when  the  source  region  is  laterally 
heterogeneous.  If  the  projection  of  first-motion  data  in¬ 
cludes  the  effect  of  propagation  through  a  plausible  velocity 
model  for  the  ocean-ridge  mantle,  the  non-orthogonality  of 
the  nodal  planes  can  be  made  to  disappear. 


k.  i 


A  particularly  clear  example  of  the  apparent  non¬ 
orthogonality  of  nodal  planes  determined  from  first-motion 
data  for  a  ridge-crest  earthquake  is  shown  in  Fig.  2.  This 
earthquake,  which  occurred  on  20  September  1969  on  the 
Reyk janes  ridge,  is  obviously  a  normal -faulting  event.  The 
angle  between  the  nodal  planes  in  the  dilatational  'quadrant' 
is  only  about  60°,  however.  Two  orthogonal  fault  planes 
could  not  be  drawn  on  the  figure  without  violating  a  large 
number  of  the  data. 

first-motion  data  in  Figure  2  were  projected  onto 
the  lower  focal  hemisphere  using  the  ray  parameter-distance 
tables  of  Herrin  (1968),  though  the  pattern  would  be  altered 
only  negligibly  if  some  other  common  earth  model  (e.g.  Jeffreys- 
Bullen)  were  used.  Of  some  effect,  though,  is  the  P-wave 
velocity  assumed  at  the  source.  We  adopted  a  value  of  8.0  km/sec 
in  plotting  this  figure.  This  value  is  reasonable  for  a  hypo- 
center  located  in  a  more  or  less  average  upper  mantle.  Assuming 
a  lower  P  velocity  appropriate  to  ocean  ridge  upper  mantle  or 
to  crustal  rocks  would  shrink  the  dilatational  first-motion 
field  even  further.  Because  we  did  not  use  such  a  lower  source- 
region  velocity,  our  estimate  of  the  nodal-plane  non-orthogon¬ 
ality  is  probably  a  conservative  one.  A  higher  P-wave  velocity 
at  the  source  would  increase  the  angle  between  the  nodal  planes 
in  Fig.  2,  but  the  velocity  would  have  to  exceed  11  km/sec  for 
the  planes  to  appear  orthogonal. 

We  cannot  conclusively  rule  out  the  possibility  that  the 
non-orthogonality  is  real,  e.g.  that  the  earthquake  source 
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is  best  represented  by  a  double  couple  and  a  superposed 
explosive  component.  An  explosive  component  with  P-wave 
amplitude  equal  to  .3  to  .5  times  the  maximum  P-wave  ampli¬ 
tude  from  the  double  couple  would  yield  the  observed  angle 
of  60°  to  70*  between  the  best  fitting  nodal  planes.  The 
nodal  surface,  would  not  be  planes  for  such  a  composite 
source,  of  course,  but  the  distinction  is  important  only 
for  propagation  directions  in  the  general  vicinity  of  the  null 
axis  on  the  focal  sphere.  Data  are  rarely  available  for  such 
paths  from  a  normal-faulting  event.  In  particular,  a  source 
consisting  of  a  double  couple  and  a  superposed  explosive  com¬ 
ponent  with  amplitude  ubout  half  the  maximum  double-couple 
ampliutde  is  entirely  consistent  with  the  first-motion  data 
in  Figure  2. 

There  are  other  methods  available  for  testing  whether  the 
nodal  plane  non-orthogonality  is  a  source  effect.  Neither 
S  waves  nor  Love  waves  generated  by  the  earthquake  should  re¬ 
flect  the  presence  of  an  explosive  component.  S-wave  polari¬ 
zation  data  for  at  least  one  ridge-crest  earthquake  (2  June 
1965)  are  not  consistent  with  the  double-couple  source  that 
would  be  obtained  by  removing  an  explosive  component  of 
sufficient  magnitude  to  account  for  the  nodal  plane  non¬ 
orthogonality.  Further,  the  seismic  moments  obtained  for 
a  ridge-crest  earthquake  using  Love  waves  would  be  expected 
to  be  systematically  lower  than  those  obtained  from  Rayleigh 
waves,  though  the  difference  (about  a  factor  of  1.5)  may  not 
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bo  discernible.  To  date,  double-couple  source  mechanisms 
have  been  adequate  to  explain  the  radiation  patterns  of 
both  Love  and  Rayleigh  waves  from  ridge  events  (Weidner  & 

Aki  1973;  Forsyth  1973)  .  In  fact,  the  source  mechanisms 
that  best  fit  the  phase  spectra  of  Rayleigh  waves  from  two 
ridge-crest  earthquakes  (Weidner  &  Aki  1973)  are  similar 
to  the  non-orthogonal  solutions  from  P-wave  first  motions 
(Sykes  1967,  1970a)  except  that  the  nodal  planes  have  a 
more  shallow  dip  and  are  orthogonal.  Another  possible  test 
might  be  made  by  examining  the  P-wave  amplitude  radiation 
pattern:  the  maximum  amplitude  of  the  compressional  lobe 

would  be  3  times  that  of  the  dilatational  lobe  for  the  sug¬ 
gested  composite  source,  though  the  critical  comparison  for 
a  normal-faulting  event  would  likely  be  between  cross-strike 
paths  to  near  stations  and  paths  through  the  earth's  core, 
a  most  difficult  comparison.  The  use  of  amplitudes  to  distin¬ 
guish  among  general  multipolar  representations  for  earthquake 
sources  has  been  discussed  by  Randall  (1972) .  For  the  reasons 
mentioned  above  and  because  of  the  clear  evidence  for  faulting 
in  the  rugged,  blocky  topography  of  the  mid-Atlantic  ridge 
(Heezen  et  al^.  1958)  we  prefer  a  simple  dislocation  model  for 
the  earthquake  mechanism  and  we  attribute  the  observed  non¬ 
orthogonality  of  nodal  planes  to  a  path  effect.  We  admit,  though, 
that  alternative  explanations  involving  the  source  cannot  yet 
be  completely  excluded. 


I 
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3 .  Velocity  Models  for  the  Ridge 

The  standard  projections  used  to  map  P— wave  observations 
onto  the  focal  sphere  are  baned  on  the  assumption  of  spherical 
symmetry.  This  assumption  is  probably  a  serious  oversimpli¬ 
fication  for  active  ridges  (and  subduction  zones)  where  rays 
may  be  strongly  bent  by  large  horizontal  gradients  of  the 
seismic  velocity. 

We  therefore  seek  a  two-dimensional  model  for  compressional— 
wave  velocity  in  the  ocean-ridge  mantle  that  can  account  for 
the  anomalous  non-orthogonality  of  nodal  planes  for  ridge-crest 
earthquakes.  This  is  in  no  sense  a  formal  inversion;  we  recog¬ 
nize  the  profound  nonuniqueness  of  the  problem.  Nonetheless, 
we  can  identify  sevei al  features  likely  to  be  of  importance 
for  a  successful  model:  ( i )  pronounced  lateral  temperature 

gradients  and  (ii)  non-horizontal  phase  boundaries,  particularly 
those  associated  with  melting.  Our  approach  is  to  construct 
models,  insofar  as  possible,  from  first  principles.  We  select 
possible  temperature  models  for  the  ocean-ridge  mantle  from 
the  literature,  and  combine  these  with  plausible  phase  diagrams 
for  mantle  material  and  measured  or  estimated  values  for  the 
dependence  of  seismic  velocity  on  temperature  and  on  phase 
(cf.  Forsyth  &  Press  1971). 

For  the  temperature  field  beneath  the  ridge  we  have  used 
the  spreading  slab  models  of  Sleep  (1974)  and  the  fluid 
dynamical  models  of  Andrews  (1972),  though  there  are  a  number 
of  quite  similar  models  in  the  literature.  Both  types  of 
models  have  some  flexibility.  For  spreading  slab  models,  we 
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adopted  values  for  slab  thickness,  temperature  at  the  slab 

base,  and  thermal  diffusivity  of  100  km,  1200°C,  and  6.6  x  10“3 
2 

cm  /sec,  respectively,  after  Sleep  (1974) .  The  only  remaining 
free  parameter  is  then  the  spreading  rate.  Andrews  (1972) 
computed  a  convection  model  for  a  spreading  half-rate  appro¬ 
priate  to  the  north  Atlantic  (1.2  cm/yr) ,  though  following 
his  suggestion  (D.J.  Andrews,  personal  communication,  1972) 
we  consider  his  temperature  field  arbitrary  to  within  multi¬ 
plication  by  a  constant  factor  not  too  different  from  unity. 
(This  is  because  the  temperature  field  scales  approximately 
as  the  thermal  conductivity  assumed  for  the  oceanic  lithos¬ 
phere)  . 

Phase  boundaries  in  the  ocean-ridge  mantle  we  estimate 
by  combining  a  chosen  temperature  field  with  Wyllxs's  (1971) 
phase  diagram  for  peridotite  in  the  presence  of  traces  of 
water  (see  Fig.  3).  Particularly  important  for  the  discussion 
below  is  the  distinction  (Ringwood  1969;  Wyllie  1971)  between 
two  regimes  of  partial  melting;  (i)  incipient  melting,  where 
temperatures  lie  between  the  wet  solidus  and  dry  solidus  for 
mantle  material  and  where  only  a  small,  temperature-insensitive 
volume  fraction  of  melt  is  present  (with  the  melt  fraction 
governed  primarily  by  water  content);  and  (ii)  normal  melting, 
at  temperatures  in  excess  of  the  anhydrous  solidus,  where 
the  melt  fraction  is  large  and  temperature  dependent. 

A  velocity  model  may  be  constructed  from  the  temperature- 
phase  field  by  adopting  values  for  the  temperature  dependence 
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of  velocity  within  single-phase  regions  and  for  the  change 
in  velocity  across  phase  boundaries.  Adequate  laboratory 
measurements  are  available  for  the  former  quantity;  we 
generally  used  a  value  appropriate  to  magnesium-rich  olivine, 
dVp/dT  =  -4  x  10  *  km/sec/C0  (Kumazawa  &  Anderson  1969)  . 

For  the  latter,  a  few  pertinent  measurements  (Spetzler  & 
Anderson  1968)  and  the  observation  that  melting  ooundaries 
are  associated  with  large  changes  in  the  shear  modulus  and 
much  smaller  changes  in  the  bulk  modulus  and  in  the  density 
allow  rough  estimates  to  be  made. 

Aside  from  the  requirement  that  the  seismic  velocity 
field  beneath  the  ridge  explain  the  several  mentioned  anoma¬ 
lous  characteristics  of  focal  mechanisms  of  ridge-crest 
earthquakes,  an  additional  strong  constrair 1  on  the  tempera¬ 
ture-phase  field  and  thus  on  the  velocity  structure  near 
the  ridge  is  provided  by  recent  observations  of  shear-wave 
attenuation  in  the  mantle  below  the  mid-Atlantic  ridge 
(Solomon  1973) .  The  variation  of  attenuation  with  propaga¬ 
tion  direction  from  the  focus  of  an  earthquake  on  the  Gibbs 
fracture  zone  and  from  the  ridge-crest  event  of  2  June  1965 
indicated  the  presence  of  a  highly  attenuating  (Q  -  10  or 
less)  zone  beneath  the  axis  of  the  ridge.  The  zone,  between 
50  and  100  km  in  width  and  no  deeper  than  50  to  150  km,  may 
be  readily  identified  as  a  region  of  extensive  melting  at 
temperatures  in  excess  of  '.he  anhydrous  solidus  of  mantle 
material  (Solomon  1973) . 
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In  Figure  4  are  shown  the  temperature  field  and  associa¬ 
ted  phase  boundaries  which  we  used  for  constructing  one 
such  velocity  mold  for  the  mid-Atlantic  ridge.  The  lateral 
ertical  extent  of  the  region  in  which  temperatures  are 
greater  than  Wyllie's  (1971)  dry  peridotite  solidus  (40  km 
haif-width,  5  to  35  km  depth)  are  in  adequate  agreement  with 
the  dimensions  of  the  region  of  very  low  Q  inferred  by  Solomon 

(1973) .  The  velocity  model  computed  from  this  temperature 
field  is  also  shown . 


This  velocity  model  was  determined  by  use  of  a  physically 
reasonable  but  admittedly  arbitrary  conversion  from  tempera¬ 
ture  and  phase  to  compressional  wave  velocity  (Fig.  3).  The 
P„  velocity  was  taken  to  be  8.15  km/sec,  the  mean  upper  mantle 
velocity  in  the  Pacific  and  Indian  Oceans  (Shor  4  Raitt  1969) . 

At  the  wet  solidus,  a  sharp  velocity  drop  (Spetrler  and  Anderson 
1968)  oi  25  km/sec  was  assumed.  Clearly  this  velocity  drop 
will  depend  on  precisely  how  much  melt  (i.e.  how  much  water) 
is  present  and  now  the  melt  is  distributed  (Walsh  1969) .  The 
quantity  of  melt  likely  at  temperatures  tetween  the  wet  and 
■Iry  solidus  is  on  the  order  of  1  percent  by  volume  (Wyllie 
1971) .  The  identification  of  the  normal  oceanic  asthenosphere 
as  a  similar  region  of  incipient  melting  (Kushiro  et  al.  1968; 
Lambert  6  Wyllie  1968;  Ringwood  1969;  Solomon  1972,  1973)  thus 
provides  the  best  guide  for  the  P-wave  velocity  at  such  temper¬ 
atures,  perhaps  7.7  to  7.9  km/sec.  Following  Wyllie  (1971), 
we  take  the  melt  concentration  to  be  approximately  independent 
of  temperature  between  the  wet  and  dry  solidus. 
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The  melt  fraction  is  a  rapidly  increasing  function  of 
temperature  above  the  dry  solidus.  We  assumed  a  smooth 
drop  in  velocity  with  increasing  temperature  in  this  melting 
range.  This  smooth  decrease  reflects  both  our  guess  of  the 
increase  in  melt  content  with  temperature,  based  on  a  compro¬ 
mise  between  estimates  of  Green  (1971)  and  Wyllie  (1971) , 
and  our  need  to  avoid  large  discontinuities  in  velocity  in  the 
ray-tracing  computations.  The  original  temperature  models  we 
used,  too,  did  not  take  account  of  the  heat  of  fusion,  so  that 
the  melt  fraction  appropriate  to  the  thermal  model  of  Fig .  4 
should  vary  smoothly  beneath  the  ridge  even  if  melt  content 
increases  sharply  with  temperature  at  the  dry  solidus  (Wyllie 
1971).  At  high  melt  concentrations,  the  S-wave  velocty  should 
be  a  function  only  of  the  bulk  modulus  and  density  of  the  solid 
and  liquid  phases.  (The  definition  of  'high  melt  concentration* 
is  arbitrary;  we  let  the  drop  in  velocity  above  the  dry  solidus 
follow  a  cosine  curve,  with  a  half  cycle  spread  over  150°C, 
in  constructing  the  model  in  Fig.  4.)  We  generally  assumed 
that  neither  bulk  modulus  nor  density  changes  significantly 
upon  melting  in  the  mantle,  so  that  the  P-wave  velocity  should 
not  drop  much  below  6  km/sec  even  in  a  region  of  extensive 
melting.  This  value  may  be  somewhat  high.  Laboratory  measure¬ 
ments  at  zero  pressure  of  compressional  wave  velocity  in 
molten  volcanic  rocks  yield  values  of  2  to  3  km/sec  (Murase 
&  Suzuki  1966;  Murase  &  KcBirney  1973) ,  though  the  implied 


4  to  1  ratio  of  compressibility  of  liquid  to  that  of  solid 
cannot  persist  to  very  high  pressure.  At  pressures  in 
excess  of  5  to  10  kbar,  the  velocity  of  the  molten  rocks  is  pro¬ 
bably  substantially  greater  than  3  km/sec. 

The  P-wave  velocity  model  in  Fig  4  is  oversimpli¬ 
fied.  No  crust  is  included,  nor  are  solid-solid  phase 
changes  (breakdown  of  plagioclase  to  garnet)  at  subsolidus 
temperatures.  Neither  of  these  omissions  are  serious, 
however.  At  most  they  introduce  horizontal  or  nearly  hori¬ 
zontal  discontinuities  in  velocity  which  will  not  affect 
the  calculations  below  or  the  conclusion  that  pronounced 
lateral  velocity  gradients  are  required  in  the  ocean-ridge 
mantle. 

4 .  Ray  Tracing  Calculations 

To  investigate  the  effect  on  seismic  wave  propagation 
of  such  two-dimensional  velocity  models  as  that  of  Fig.  4, 
we  resort  to  geometric  ray  theory.  While  ray  theory  may 
have  questionable  validity  for  the  long  periods  that  dominate 
the  body-wave  spectra  from  ridge-crest  earthquakes,  it  is 
probably  adequate  for  the  higher  frequencies  which  control 
P-wave  first  motions.  The  computational  scheme  by  which  rays 
are  traced  in  three  dimensions  through  models  such  as  in 
Fig.  4  has  been  fully  documented  by  Julian  (1970)  and  has 
been  applied  to  laterally  heterogeneous  models  of  subduction 
zones  by  Toksflz,  Minear  and  Julian  (1971)  and  by  Davies  and 
Julian  (1972) . 
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Ray  paths  from  a  surface  source  on  the  ridge  axis  of 
the  velocity  model  of  Fig.  4  are  shown  in  Fig.  5.  All  of 
these  rays  show  a  pronounced  downward  curvature.  Rays 
leaving  the  source  with  take-off  angles  (measured  from 
vertical)  less  than  15°  are  bent  nearly  to  the  vertical  or 
beyond.  Rays  leaving  the  source  at  take-off  angles  of  40° 
to  50°  are  all  bent  downward  by  some  16°.  In  effect,  the 
pocket  of  low-velocity  material  associated  with  extensive 
melting  acts  as  a  cylindrical  lens,  focusing  most  downward 
traveling  waves  toward  the  vertical. 

In  general,  both  the  horizontal  phase  velocity  and 
azimuth  of  a  ray  are  affected  by  propagation  through  the 
laterally  heteroegeneous  velocity  structure.  A  convenient 
graphical  description  of  these  effects  is  shown  in  Fig.  6. 
This  figure  may  be  regarded  as  a  correction  to  the  standard 
projections  used  to  map  teleseismic  P-wave  observations 
back  tc  the  focal  sphere.  The  tail  of  each  arrow  in  Fig.  6 
indicates  the  take-off  angle  (or,  alternatively,  the  ray 
parameter)  and  azimuth  of  a  ray  leaving  the  source.  The  head 
of  the  arrow  is  the  location  on  the  focal  sphere  of  the  same 
ray  that  we  would  have  naively  predicted  had  we  measured 
the  ray  parameter  and  azimuth  of  that  ray  after  propagation 
through  the  velocity  model.  The  propagation  direction 
represented  by  the  head  of  the  arrow  governs  the  epicentral 
distance  at  which  the  ray  will  airive  at  the  earth's  surface 
and  the  apparent  velocity  that  a  seismic  array  would 
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measure.  It  is  the  propagation  direction  represented  by  the 
tail  of  the  arrow  that  is  appropriate,  however,  for  plotting 
first  motions  on  the  focal  sphere. 

The  overall  effect  of  the  ridge-mantle  structure  on 
body-wave  propagation  directions  is,  according  to  Fig.  6, 
quite  pronounced.  All  rays  are  bent  downward  and  toward  the 
ridge  axis,  with  the  amount  of  bending  very  sensitive  to 
take-off  angle  and  azimuth.  The  model  predicts  that  some 
waves  wi.1 1  cross  beneath  the  ridge  axis,  producing  a  region 
of  multiple  arrivals  (triplication);  this  phenomenon  is 
particularly  evident  for  rays  leaving  the  source  at  azimuths 
nearly  along  the  ridge  and  at  high  incidence  angles. 

The  P-wave  first  motions  from  the  ridge-crest  earthquakes 
of  20  September  1969,  with  the  position  of  each  point  on  the 
focal  sphere  corrected  according  to  Fig.  6,  are  shown  in 
Fig.  2.  With  the  corrections,  all  points  are  moved  outward 
from  the  ridge  axis  and  two  orthogonal  nodal  planes  can  be 
fit  to  the  data  with  no  difficulty.  (The  corrections  of 
Fig.  6,  although  generally  appropriate  to  much  of  the  north 
Atlantic,  are  not  strictly  applicable  to  the  obliquely  spread¬ 
ing  Reyk janes  ridge.  A  slightly  slower  spreading  rate  would 
not  affect  the  corrections  appreciably,  however;  see  Fig.  9 
below. ) 

It  should  be  noted  that  even  if  the  velocity  model  of 
Fig.  4  were  a  correct  representation  of  the  mantle  beneath 
the  mid-Atlantic  ridge,  Fig.  6  is  appropriate  only  for  a 
particular  epicenter  and  focal  depth.  There  is  some  evidence 
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that  ridge  crest  earthquakes  may  be  shallower  than  a  few  km 
in  depth  (Weidner  &  Aki  1973) .  The  present  uncertainty  in 
epicentral  coordinates  of  ridge  crest  earthquakes,  however, 

is  about  10  to  20  km,  and  normal  faulting  earthquakes  further 
than  50  km  from  the  ridge  crest  have  been  observed  (Sykes  1970b) . 

We  therefore  investigated  the  effect  on  ray  paths  of  moving 
the  surface  source  off  the  ridgn  crest  (Fig.  7).  For  sources 
within  several  tens  of  kilometers  of  the  ridge  axis,  the  effect 
on  propagation  directions  of  the  velocity  structure  beneath  the 
ridge  remains  pronounced,  though  strongly  dependent  on  the 
source-axis  offset.  Both  triplications  and  modest  shadow  zones 
are  evident  in  the  ray  tracings  of  Fig.  7  for  sources  located 
5,  10  and  20  km  from  the  ridge  axis.  For  sources  more  distant 
from  the  axis  (e.g.  100  km),  che  effect  on  the  velocity 
structure  is  very  modest. 

5 .  Discussion  of  Velocity  Models 

The  velocity  structure  of  Fig .  4  is  just  one  velocity 
model,  based  on  a  reasonable  temperature  distribution,  that 
can  provide  an  explanation  for  the  apparent  non-orthogonality 
of  nodal  planes  for  ridge  crest  earthquakes.  The  critical 
question  is  which  aspects  of  the  model  are  necessary  for  such 
an  explanation  to  hold  and  which  aspects  may  be  varied. 

We  investigated  a  large  number  of  similarly  constructed 
velocity  models,  with  different  temperature  models  and  differ¬ 
ent  rules  for  converting  from  temperature  and  phase  to  seismic 
velocity.  The  universal  characteristic  of  all  models  that  success¬ 
fully  explained  the  non-crthogonality  of  Fig.  2  was  a  pro¬ 
nounced  lateral  velocity  gradient;  changes  in  velocity  of 
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1  km/sec  over  horizontal  distances  of  15  to  25  km  were  gener¬ 
ally  required  at  some  fairly  shallow  depth  in  the  velocity 
model.  A  lateral  change  of  such  a  magnitude  cannot  arise  from 
temperature  variations  in  material  that  is  everywhere  solid 
and  mineralogically  homogeneous.  Similarly,  such  a  variation 
almost  certainly  cannot  be  due  to  a  solid-solid  phase  change 
in  the  lithosphere,  both  became  the  velocity  discontinuities 
at  phase  changes  likely  at  such  depths  are  too  small  and  be¬ 
cause  the  pressure  derivative  of  the  transition  temperature 
for  such  reactions  are  generally  too  large.  Rather  such  a 
horizontal  velocity  gradient  must  be  associated  with  melting 
beneath  the  ridge,  melting  of  such  an  extent  to  produce  P-wave 
velocities  near  7  km/sec  or  less  in  mantle  material. 

There  is  some  arbitrariness , however , in  the  parameters  that 
specify  the  shape  of  the  molten  region,  particularly  the  depth 
to  the  top  of  the  region  of  extensive  melting  and  the  depth  in¬ 
terval  over  which  such  extensive  melting  is  prevalent. 

As  an  example ,  in  Fig.  8  are  shown  the  ray  paths 
from  a  surface  source  in  a  P-wave  velocity  model  foi.  the  ocean- 
ridge  mantle  based  on  the  temperature  calculations  of  Andrews 
(1972).  Specifically  the  temperature  distribution  is  Andrews' 
model  for  the  north  Atlantic,  scaled  by  a  factor  cf  0.54  (see 
above  and  Solomon  1973)  .  Temperatures  in  the  upper  20  km  be¬ 
neath  the  ridge  are  cooler  in  this  model  than  in  that  of  Fig.  4, 
but  temperatures  at  100  km  depth  are  about  100°C  hotter.  The 
simple  phase  diagram  of  Fig.  3  then  predicts  extensive  melting 
in  a  region  15  to  50  km  deep  extending  to  a  distance  of  roughly 
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60  Ion  from  the  ridge  axis .  a  velocity  model  was  generated  from 
the  temperature-phase  field  in  the  manner  of  the  model  of  Fig.  4, 
except  that  the  decrease  in  P-wave  velocity  with  temperature 
above  the  anhydrous  solidus  was  taken  to  occur  over  a  narrower 
temperature  interval  (60°C) . 

Such  a  velocity  model  also  yields  the  requisite  downward 
bending  of  rays  (Fig.  8) .  The  dependence  on  take-off  angle  is 
somewhat  different  than  in  Fig.  5, however .Rays  leaving  the  source 
at  take-off  angles  less  than  15*  are  bent  only  slightly.  The 
maximum  effect. is  for  rays  with  take-off  angles  in  the  range  40* 
to  60°;  these  rays  are  bent  downward  by  some  15°  to  20°.  If 
the  appropriate  corrections  to  focal  sphere  positions  based  on 
this  model  were  applied  to  the  first  motions  of  Fig.  2,  a  source 

mechanism  with  orthogonal  nodal  planes  could  be  readily  fit  to 
the  data. 

Thus  the  ray  tracing  calculations  provide  no  firm  constraint 
the  depth  to  the  top  of  the  region  of  extensive  melting.  Tem¬ 
peratures  exceed  the  anhydrous  solidus  at  depths  greater  than 
5  km  in  the  model  of  Fig.  4;  they  exceed  the  anhydrous  solidus 
at  depths  greater  than  15  km  in  the  model  on  which  Fig.  8  was 
based;  and  both  models  are  acceptable  from  the  standpoint  of 
explaining  the  anomalous  non-orthogonality  of  nodal  planes  for 
ridge-crest  earthquakes.  At  present,  other  types  of  observations 
offer  the  best  hope  for  discriminating  among  seismic  models. 

For  instance,  magnetotelluric  data  suggest  that  th»  anhydrous 
solidus  of  peridotite  is  20  km  deep  beneath 
Iceland  (Hermance  &  Grillot  1970).  Also,  the  basaltic  rocks 
of  the  ocean  floor  are  thought  to  be  the  product  of  extensive 
partial  melting  and  rapid  separation  of  liquid  from  solid 
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at  15  to  25  km  depth  beneath  mid-ocean  ridges  (Kay,  Hubbard  & 
Gast  1970) .  We  should  note  in  this  context  that  the  low  seismic 
velocities  in  our  models  cannot  per  se  distinguish  extensive 
in  situ  partial  melting  from  an  aggregation  of  melt  that  has 
separated  from  deeper  regions  of  lesser  melt  content.  Geoche¬ 
mical  evidence  that  ocean-floor  basalts  are  the  products  of 
extensive  (10  to  30  percent  by  volume)  partial  melting  (Kay  et 
al.  1970)  and  Sleep's  (1974)  arguments  on  the  inefficiency  of 
magma  segregation  at  low  melt  concentrations  convince  us,  how¬ 
ever,  that  our  temperature-phase-velocity  models  are  not  un¬ 
reasonable.  Some  upwards  segregation  of  melt  immediately  be¬ 
neath  the  ridge  axis  (Atwater  &  Mudie  1973)  may  deplete  the 
melt  fraction  at  some  distance  from  the  axis  over  what  would 
be  estimated  from  temperature  and  the  phase  diagram  of  un¬ 
differentiated  mantle  material. 

The  model  calculations  do  provide  some  constraint  on  the 
lateral  extent  of  melting.  In  the  north  Atlantic,  at  least, 
the  region  of  high  melt  content  (temperatures  above  the  anhy¬ 
drous  solidus)  cannot  be  confined  to  a  narrow  dike  (cf.  Kay 
et  al.  1970) ,  nor  does  it  likely  extend  to  several  hundred 
kilometers  from  the  ridge  axis  (cf .  Oxburgh  &  Turcotte  1968) . 
For  lithosphere  spreading  rates  comparable  to  that  in  the 
north  Atlantic,  at  any  given  depth  in  the  range  15  to  40  km 
the  phase  boundary  identified  with  the  anhydrous  peridrtite 
solidus  extends  from  the  ridge  axis  a  distance  of  roughly  one 
to  two  times  the  depth.  The  maximum  half -width  of  the  region 
of  extensive  melting  is  some  40  to  60  km.  That  region  of 
extensive  melting  is  relatively  broad,  in  contrast 
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to  the  narrow  zone  of  active  intrusion  on  the  ridge  axis 
(Atwater  &  Mudie  1973) ,  argues  that  a  significant  quantity 
of  basaltic  material  is  >dded  to  the  cooling  oceanic  litho¬ 
sphere  below  the  crust  (Press  1969;  Forsyth  &  Press,  1971)  , 

To  summarize,  velocity  models  for  the  ocean-ridge  mantle 
that  explain  the  unusual  first-motion  patterns  of  normal- 
faulting  earthquakes  on  the  ridge  crest  all  involve  substan¬ 
tial  horizontal  velocity  gradients  which  can  only  be 
explained  by  partial  melting  at  shallow  depths  beneath 
spreading  centers.  An  obvious  corollury  is  that  the  details 
of  the  conventional  fault-plane  solution  of  any  earthquake 
on  the  mid-ocean  ridge  system  are  suspect.  The  nature  of 
faulting  (e.g.  normal  Dr  strike-slip)  will  not  be  in  doubt 
(Sykes  1967)  but  such  important  particulars  as  slip  vectors 
or  inferred  principal  axes  of  pre-earthquake  stress  may  well 
be.  This  will  be  especially  true  of  sources  with  non-vertical 
fault  planes  (e.g.  most  normal  faults),  but  may  also  be  true 
to  a  lesser  extent  for  strike-slip  earthquakes  located  on 
transform  faults  near  the  intersection  with  a  ridge  crest. 

6 .  Further  Possible  Tests 

There  are  a  number  of  means  by  which  the  general  features 
of  the  velocity  models  presented  above  can  be  tested  and 
through  which  the  range  of  possible  models  can  be  reduced. 

A  logical  extrapolation  of  the  above  discussion  is  to 
include  S-wave  polarization  data.  As  a  result  of  propagation 
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through  a  laterally  heterogeneous  velocity  model,  both  the 
propagation  direction  and  the  polarization  angle  will  change. 

If  S-wave  polarization  data  are  corrected  for  such  changes, 
they  should  then  be  consistent  with  the  fault-plane  solution 
deduced  from  corrected  P-wave  first  motions.  The  S-wave 
velocity  structure  can  be  calculated  in  a  manner  analogous 
to  that  used  above,  and  should  serve  as  a  check  on  the  assump¬ 
tions  employed  in  constructing  the  P-vclocity  model,  particu¬ 
larly  the  behavior  of  bulk  modulus  and  rigidity  at  temperatures 
in  excess  of  the  solidus.  A  potential  difficulty,  however, 
is  that  the  S-wave  polarization  angles  from  ridge-crest  earth¬ 
quakes  must  be  measured  on  long-period  waves,  for  which  geo¬ 
metric  ray  theory  may  be  inadequate. 

The  models  discussed  above  have  been  tailored  specifically 
to  ridges  in  the  north  Atlfu'tic.  A  legitimate  question  is 
whether  the  procedures  usrd  to  derive  these  models  have  general 
applicability  to  spreading  centers  throughout  the  world.  We 
therefore  examined  a  number  of  models  similar  in  every  detail 
to  that  in  Fig ,  4  except  that  the  spreading  rate  was 
varied.  Ray  paths  from  surface  sources  on  the  axes  of 

ridges  spreading  at  a  number  of  rates  are  shown  in  Fig.  9. 

The  thermal  model  calculated  according  to  Sleep  (1974) 
for  a  spreading  half-rate  of  0.5  cm/yr  is  somewhat  cooler 
beneath  the  ridge  than  faster  spreading  models.  Temper¬ 
atures  in  excess  of  the  dry  solidus  are  reached  only  within 
a  narrow  dike  immediately  beneath  the  ridge  axis.  The  lateral 
gradients  in  P-wave  velocity  associated  with  this  temperature- 
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phase  distribution  are  thus  more  modest  than  those  in  Fig.  4 
but  are  still  significant.  Some  downward  bending  of  rays 
is  evident  for  this  model  in  Fig.  9.  Rays  with  initial  take¬ 
off  angles  of  30*  to  50*  are  bent  down  6°  to  7°,  so  that  a 
pure  dip-slip  normal  fault  with  each  nodal  plane  in  reality 
dipping  at  45°  to  the  horizontal  would  appear  to  have  nodal 
planes  non-orthogonal  by  about  13°.  This  is  in  approximate 
accord  with  observed  non-orthogonality  of  nodal  planes  for 
ridge-crest  earthquakes  in  the  Arctic  (Table  1}  and  is  not 
inconsistent  with  less  restrictive  first  motion  data  from  earth¬ 
quakes  on  other  ridges  with  comparable  spreading  rates  (e.g. 
the  Red  Sea,  McKenzie,  Davies  &  Molnar  1970) . 

For  ridges  with  spreading  half-rates  between  perhaps  0.5  and 
about  2.0  cm/yr,  the  apparent  non-orthogonal ity  of  nodal  planes 
predicted  for  surface  normal^-faulting  earthquakes  on  the  ridge 
axis  increases  with  spreading  rate.  Thus  for  the  model  in 
Fig.  9  with  a  half -rate  of  1.0  cm/yr,  the  downward  bending  of 
rays  is  quite  similar  but  somewhat  less  than  that  shown  in  Fig.  5. 
If  the  first  six  earthquakes  in  Table  1  are  arranged  in  order  of 
increasing  spreading  rate,  then  the  reported  angle  between 
nodal  planes  decreases  about  as  regularly  as  one  might  expect 
considering  the  simplifications  in  the  models,  the  uncertainty 
in  the  determination  of  each  nodal  plane,  and  the  uncertainty 
in  hypocentral  coordinates  of  each  earthquake. 

The  most  pronounced  distortion  of  ray  paths  occurs  for  spread¬ 
ing  half-rates  near  2.0  cm/yr  (Fig.  9).  Rays  with  initial 
take-off  angles  of  30°  to  50°  are  bent  towards  the  downward 


verbical  by  21°  to  23°.  If  this  model  is  correct,  therefore, 
fault-plane  solutions  of  normal  faulting  earthquakes  on 
radges  with  half  rates  of  2  cm/yr  should  show  nodal  planes 
at  an  apparent  angle  of  roughly  45®.  No  such  phenomenon  has 
been  observed,  but  neither  do  the  currently  scanty  first-motion 
data  from  earthquakes  on  such  ridgos  rule  out  such  extreme 
effects.  The  angle  between  nodal  planes  listed  in  Table  1 
for  the  event  in  the  Gulf  of  California  studied  by  Thatcher 
&  Brune  (1971)  is  something  of  an  upper  bound;  a  lower  bound 
on  this  angle  would  be  considerably  smaller  than  45®.  Portions 
of  the  Carlsberg  ridge  are  spreading  at  a  half  rate  near  2.0  cm/yr. 
P-wave  first  motions  from  ridge-crest  earthquakes  in  that  region, 
though  interpreted  (Banghar  &  Sykes  1969)  in  terms  of  orthogonal 
nodal  planes,  either  provide  no  constraint  on  one  or  both  planes 
or  are  better  fit  by  planes  non-orthogonal  by  several  tens  of 
degrees  (see  especially  their  Fig.  4) . 

At  high  spreading  rates  (5  cm/yr  half -rate  or  more),  the 
zone  of  extensive  melting  is  very  broad,  so  the  implied  hori¬ 
zontal  gradients  in  P-wave  velocity  beneath  the  ridge  are  small 
except  at  the  shallow  and  nearly  horizontal  lithosphere- 
asthenosphere  boundary.  For  such  velocity  models,  the  effect 
on  ray  paths  is  less  severe  than  at  lower  spreading  rates,  and 
the  downward  bending  increases  monotonically  with  initial  take¬ 
off  angle.  Rays  with  starting  take-off  angles  in  the  range 
20°  to  50®  are  bent  downward  by  5®  to  13®  in  the  model  with 
5  cm/yr  spreading  half-rate  (Fig.  9) .  These  changes  are 
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reduced  by  a  factor  of  three  in  a  model  (not  shown)  with 
10  cm/yr  spreading  half-rate.  To  our  knowledge,  there  are 
no  normal-faulting  earthquakes  on  such  rapidly  spreading 
ridges  large  enough  to  determine  a  fault-plane  solution  from 
first  motions. 

The  more  traditional  seismological  tools  for  deciphering 
velocity  structure  may  also  serve  as  additional  tests  or 
constraints  on  models  such  as  Fig.  4.  Long  refraction  lines, 
out  to  distances  great  enough  so  that  waves  bottoming  below 
the  asthenosphere  are  first  arrivals,  would  provide  the  most 
straightforward  test.  Potential  problems  are  the  very  low  Q 
in  the  asthenosphere  beneath  the  ridge  axis  and  the  likelihood, 
if  Fig  6  can  be  a  guiui,  of  a  confusing  assortment  of  laterally 
refracted  arrivals.  Francis  (1969)  has  constructed  such  a 
long  refraction  profile  along  the  mid-Atlantic  ridge  using 
earthquake  sources  and  pairs  of  seismograph  stations  cn  Iceland. 
His  apparent  velocity  data  showed  a  good  deal  of  scatter;  appar¬ 
ent  velocities  of  P  waves  as  low  as  6  km/sec  at  7°  epicentral 
distance  were  recorded.  Nonetheless  the  models  that  fit  his 
smoothed  data  bear  some  resemblance  to  the  velocity-depth  dis¬ 
tribution  beneath  the  ridge  axis  in  Fig.  4.  A  long  marine 
refraction  profile  using  artificial  sources  and  probably  using 
bottom  seismometers  in  addition  to  hydrophones  would  be  the 
most  fruitful  way  to  test  the  details  of  the  ccean-ridge  mantle 
velocity  structure. 

Inversion  of  'pure-path'  surface-wave  velocities  also  holds 
promise  for  resolving  some  of  the  details  of  lateral  variation 
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seismic-wave  velocity  within  the  lithosphere  and  astheno- 
sphere  (Weidner  1974;  Forsyth  1973).  The  increase  of  litho¬ 
sphere  thickness,  identified  as  the  depth  to  the  wet  solidus, 
with  age  in  Fig.  4  is  in  approximate  agreement  with  Weidner' s 
(1974)  results  but  is  less  than  in  Forsyth's  (1973)  models.  The 
differences  in  Rayleigh-wave  velocities  and  inferred  models 
between  these  studies  may  be  real  or  may  be  a  result  of  severe 
lateral  refraction  effects  in  the  Atlantic  Ocean  data  (Forsyth 
1973).  The  growth  of  the  lithosphere  with  time  in  Fig.  4  could 
be  brought  into  better  agreement  with  Forsyth's  (1973)  models 
by  reducing  the  estimated  water  content  of  the  peridotite  mantle 
and  thus  raising  the  wet  solidus  (Krshiro  et  al.  1968) .  Because 
the  surface-wave  velocities  are  more  sensitive  to  shear  velocity 
than  to  compressional  velocity  and  because  surface 
waves  average  laterally  over  distances  comparable  to  their  wave¬ 
length,  the  details  of  models  such  as  Fig.  4,  especially  the 
region  of  extensive  melting  and  very  low  velocities  beneath  the 
ridge  crest,  are  not  well  suited  to  study  by  surface-wave 
analysis. 

A  simple  measure  of  an  anomalous  velocity  structure  such 
as  in  Fig.  4  is  the  travel-time  delay  of  seismic  waves  which  have 
propagated  through  such  a  region.  If  one  knew  the  origin 
time  and  hypocenter  of  a  ridge-crest  event,  then  measurement 
of  the  delay  (with  respect  to  standard  tables)  as  a  function 
of  propagation  direction  would  be  straightforward.  Thatcher 
&  Brunt!  (1971)  have  estimated  in  this  manner  the  travel-time 
delay  af  P-waves  from  events  on  the  spreading  center  in  the 


Gulf  of  California  to  be  about  2  sec.  Conversely,  the  average 
delay  at  a  station  on  the  ridge  crest  could  be  compared  to 
that  at  a  station  overlying  older  sea  floor  if  records  from 
a  number  of  teleseismic  sources  were  available.  The  P-wave 
delay  at  Iceland,  for  instance,  is  over  1  sec  greater  than  at 
stations  in  stable  continental  areas  of  Greenland,  Scandinavia, 
and  Scotland  (Tryggvason  1964;  Long  &  Mitchell  1970). 

As  a  guide  to  possible  future  studies  of  delay  times  at 
rid^c  crests,  we  calculated  the  P— wave  delays  for  a  ridge- 
crest  source  atop  the  velocity  model  of  Fig.  4  relative  to 
the  * t avel— times  for  a  source  at  the  edge  of  our  velocity  model 
(roughly  20  m.y.  old  sea-floor) .  These  delays  are  shown 
plotted  on  a  focal-sphere  projection  (Davies  &  McKenzie  1969) 
in  Fig.  10.  The  mean  residual  averaged  over  the  focal  sphere 
is  about  +0.6  sec.  This  delay  would  be  slightly  greater  if  the 
comparison  were  made  with  lithosphere  older  than  20  m.y.  The 
variation  of  residual  with  propagation  direction  is  slight, 
except  for  paths  at  large  take-off  angles  nearly  along  the 
ridge  axis.  In  this  situation,  there  is  an  ambiguity  between 
the  travel  time  anomaly  and  the  unknown  origin  time  of  an 
earthquake.  Thus,  teleseismic  travel  times  could  be  used  to 
constrain  the  velocity  structure  of  the  ocean-ridge  mantle 
only  if  the  hypocenter  and  origin  time  of  a  ridge-crest  event 
were  known  from  independent  information;  and  even  then  the 
only  constraint  would  probably  be  a  travel-time  offset  largely 
insensitive  to  propagation  direction. 


-136- 


In  conclusion/  there  are  a  number  of  additional 
techniques  which,  with  various  likelihoods  of  success,  offer 
the  potential  for  testing  the  velocity  models  proposed  to 
explain  the  apparent  non-orthogonality  of  nodal  planes  for 
ridge-crest  earthquakes.  The  most  promising  are  long  refrac¬ 
tion  profiles,  examination  of  S-wave  polarization  data  from 
similar  earthquakes,  and  possibly  the  measurement  of  travel¬ 
time  delays  of  P  and  S  waves  from  sources  on  or  near  ridge- 
crests. 
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plane  constrained  primarily  by  proximity  of  low-amplitude  first  arrivals. 


Figure  Captions 

Fig.  1  Location  of  ridge-crest  earthquakes  for  which  P-wave 
first  motions  indicate  non-orthogonal  nodal  planes. 
Events  ar .*  numbered  as  in  Table  1.  Plate  boundaries  are 
shown  by  heavy  lines.  Azimuthal  equidisvant  projection 
about  the  north  pole. 


Fig.  2  Left:  Conventional  fault-plane  solution  from  P-wave 
first  motions  for  the  earthquake  of  September  20,  1969 
(equal  area  projection) .  Open  circles  are  dilatations, 
closed  circles  are  compressions,  and  circles  with  superposed 
crosses  indicate  probable  proximity  to  a  nodal  plane. 
Smaller  symbols  represent  somewhat  less  reliable  readings 
than  larger  ones.  All  data  were  read  from  long-period 
records  of  the  WWSSN  or  the  Canadian  network.  The  strike 
$  and  dip  6  of  the  implied  nodal  planes  are  also  given. 

The  angle  between  the  nodal  planes  is  60°. 

Right:  Same  first-motion  data  corrected  for  propagation 
through  the  ocean-ridge  mantle  (Fig.  4)  according  to  Fig.  6. 
The  new  fault-plane  solution  has  orthogonal  nodal  planes. 

Fig.  3  Melting  relations  for  peridotite  in  the  presence  of 

0.1%  water,  from  Wyllie  (1971) .  The  solidus  at  low  pres¬ 
sures  (shallow  depths)  is  controlled  by  the  breakdown  of 
hornblende. 

Inset:  Adopted  dependence  of  comp.,  ess  ional  wave  velocity 

V  on  temperature  and  phase.  The  particular  curve  shown 
P 

corresponds  to  20  km  depth  but  the  sharp  drop  at  the  wet 
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Fig. 


Fig. 


Fig. 


solidus  (T^  and  the  smooth  decrease  above  the 
anhydrous  solidus  (T2)  are  maintained  at  all  depths. 

4  Top:  Isotherms  and  regions  of  melting  in  the 
spreading  oceanic  lithosphere.  The  temperature  distri¬ 
bution  is  calculated  according  to  Sleep  (1974)  for  a 
spreading  half-rate  of  1.2  cm/yr.  The  regions  of 
incipient  and  extensive  melting  are  determined  from  this 
temperature  distribution  and  the  phase  diagram  of  Fig.  3. 

Bottom:  Contours  of  constant  compressional-wave 
velocity.  The  conversion  from  temperature  and  phase  to 
velocity  follows  Fig.  3  and  the  arguments  in  the  text. 

5  Ray  paths  from  a  surface  earthquake  on  the  axis  of 
the  ridge  model  of  Fig .  4 .  All  rays  remain  in  a  vertical 
plane  perpendicular  to  the  ridge  axis.  Distances  along 
the  earth's  surface  are  marked  in  km;  no  vertical  exaggera¬ 
tion.  Rays  are  spaced  by  5*  in  initial  take-off  angle, 
from  0*  to  ±70°. 

6  Theoretical  distortion  of  apparent  locations  of  P  waves 
on  the  focal  sphere  caused  by  propagation  through  the 
upper  mantle  model  of  Fig.  4  from  a  surface  focus  on  the 
axis  of  a  north-south  ridge.  Each  arrow  connects  the 
location  of  the  ray  at  the  focus  (tail  of  arrow)  with  the 
location  corresponding  to  the  ray  direction  at  a  depth 

of  80  km  (head  of  arrow  -  indicated  by  asterisk) .  The 
lower  focal  hemisphere  is  shown  in  equal  area  projection 
and  rays  are  spaced  by  10°  in  azimuth  and  take-off  angle 
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(from  20°  to  M#)  at  the  focus.  Only  rays  in  the 
right  half  of  the  focal  sphere  are  included  for  clarity. 

Fig.  7  Hay  paths  from  surface  earthquakes  at  various  distances 
(5,  10,  20,  100  km)  from  the  axis  of  the  ridge  model  of 
Fig.  4.  Other  particulars  as  in  Fig.  5. 

Fig.  8  Ray  paths  from  a  surface  earthquake  on  the  axis  of  an 
alternative  ridge  model  based  on  the  temperature  model  of 
Andrews  (1972);  see -text.  Other  particulars  as  in  Fig.  5. 

Fig.  9  Dependence  of  distortion  of  ray  paths  on  spreading 
rate.  All  ray  paths  are  from  surface  events  on  ridge 
axes.  Velocity  models  are  calculated  analogously  to  that 
in  Fig.  4  except  for  the  different  spreading  rates.  Other 
particulars  as  in  Fig.  5. 

Fig.  10  Left:  P-wave  travel- time  residuals  predicted  for  a 

surface  source  on  the  ridge  axis  of  the  model  of  Fig.  4. 
Plusscs  are  positive  residuals;  circles  are  negative 
residuals.  The  diameter  of  a  symbol  is  linearly  propor¬ 
tional  to  the  magnitude  of  the  residual;  the  symbols  at 
lower  center  are  ±1.0  sec.  Points  are  plotted  on  the 
focal  sphere  (equal  area  projection)  according  to  the 
azimuth  and  ray  parameter  of  the  wave  at  the  bottom  of 
the  laterally  heterogeneous  ridge  model.  Only  rays 
leaving  the  source  in  the  right  half  of  the  lower  focal 

sphere  are  plotted  for  clarity. 

Right:  Same  data  after  subtracting  0.6  sec  from  each 
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5.  ARRAYS  AND  DATA  PROCESSING 

K .  1  Automatic  Phase  Identification  With  One  and  Two  Large 
Aperture  Seismic  Arrays  by  S.  Shlien  and  M.N.  ToksOz 
(abstract) 

Schemes  to  automatically  identify  later  phases  using 
LASA  and  using  LASA  together  with  NORSAR  were  designed  and 
tested.  With  a  single  array,  phase  identification  was 
accomplished  by  searching  for  a  P  or  PKP  phase  that  could 
be  associated  with  a  particular  later  phase  such  as  PP, 

PcP,  etc.  On  the  basis  of  the  time  interval  between  the 
signals,  their  slowness,  and  relative  amplitudes  the  best 
phase  interpretation  was  chosen  by  means  of  the  Neyman- 
Pearson  likelihood  ratio  test.  With  two  arrays,  detections 
from  one  array  were  checked  against  detections  from  the 
other  array  for  signals  originating  from  the  same  event. 

Testing  the  schemes  on  the  LASA  and  NORSAR  Detection 
Logs  9  later  phases  per  day  were  detected  and  11  events 
per  day  were  found  common  to  the  LASA  and  NORSAR  detection 
logs.  Confusion  of  later  phases  was  very  minor.  The  most 
serious  problem  wan  the  identification  of  fictitious  phases 
and  events  that  occurred  due  to  various  coincidences. 
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5 • 2  An  Event  Recording  Syr  tern  for  Monitoring  Small 

Earthquakes  by  Bruce  P.  Ambuter  and  Sean  C.  Solomon 
Summary 

To  monitor  small  earthquakes  on  the  ocean  bottom  or 
in  remote  areas  on  land,  a  seismic  recording  system  should 
operate  «t  low  power  for  extended  periods  of  time,  should 
reproduce  three  components  of  ground  motion  with  a  large 
dynamic  range,  and  should  preferably  be  small  and  inexpen¬ 
sive.  We  describe  such  a  system  based  on  an  event  detector, 
which  continuously  surveys  the  background  noise,  sets  a 
threshold,  and  triggers  an  inexpensive  digital  tape  recorder 
when  a  seismic  signal  exceeds  that  threshold.  Event  recording 
offers  the  advantages  over  conventional  continuous  recording 
of  more  efficient  use  of  storage  capacity  and  easier  later 
interpretation.  A  key  element  of  the  system  is  a  continuously 
updated  semiconductor  buffer  memory,  which  assures  that  first 
arrivals  are  included.  Average  power  consumption  is  about  a 
watt.  The  system  design  is  e^^ily  adaptable  to  monitoring 
other  types  of  rare  events  of  unpredictable  occurrence. 
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INTRODUCTION 

The  operation  of  a  seismometer  system  in  a  remote  or 
nearly  inaccessible  area  has  always  posed  logistical,  mechani¬ 
cs  and  financial  problems.  Field  instruments  or  temporary 
stations  must  operate  efficiently  at  low  power,  particularly 
when  continuous  monitoring  over  long  time  intervals  is  desired. 
At  the  same  time,  such  standard  seismological  tools  as  spectral 
analysis  require  that  the  data  recording  system  have  a  broad 
or  flexible  frequency  response  and  a  wide  dynamic  range.  Moti¬ 
vated  primarily  by  a  desire  to  monitor  small  earthquakes  on  the 
ocean  bottom  or  at  temporary  land  stations  with  a  movable  array 
of  three-component  seismometers,  we  developed  a  seismic  recording 
package  that  meets  these  diverse  requirements.  The  heart  of  the 
system,  described  in  this  paper,  is  an  event  recording  scheme 
based  on  a  continuously  updated  semiconductor  memory. 

There  are  two  distinctly  different  methods  for  recording 
seismic  data.  The  first  is  continuous  sto-ing  of  seismometer 
output.  The  second  is  event  recording,  i.e.,  editing  data  and 
storing  only  after  a  threshold  has  been  exceeded.  The  first 
system  is  more  simple  to  design.  The  second  requires  less 
storage  capacity  and  makes  later  reading  and  interpreting  of 
the  data  loss  time-consuming. 

Conventional  seismometer  systems  with  adequate  dynamic 
range  for  spectral  analysis  generally  use  continuous  recording, 
either  in  the  FM  (frequency-modulation)  or  direct  analog  mode. 
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Fach  mode  has  both  advantages  and  drawbacks.  Direct  analog 
recording  is  the  easiest  system  to  implement,  requiring  a 
minimum  of  components.  Its  main  disadvantages  are  susceptibi¬ 
lity  to  tape  dropouts  and  non-linear  low  frequency  response. 

With  FM  recording,  there  is  a  more  uniform  frequency  response. 
However,  there  is  a  trade-off  between  tape  speed  and  high  fre¬ 
quency  response  that  can  unduly  limit  recording  time.  In 
addition,  the  recording  speed  stability  is  a  critical  parameter. 
In  either  system,  the  dynamic  range  is  limited  to  about  40  db 
and  the  tape  recorder  is  commonly  large  and  expensive.  None 
of  these  disadvantages  are  fatal,  we  should  note,  and  many  can 
be  circumvented  by  careful  design  work  (Dibble,  1964;  Sacks, 
1966;  Muirhead  and  Simpson,  1972;  Green,  1973).  Continuous 
recording  can  also  be  performed  with  a  digital  recorder,  though 
the  size  and  cost  of  such  a  unit  must  necessarily  be  large 
for  most  recording  times  of  interest. 

An  event  recorder  which  continuously  surveys  the  ambient 
noise  at  a  recording  site,  sets  a  threshold,  and  triggers  the 
recorder  when  a  seismic  signal  exceeds  the  threshold,  requires 
far  less  storage  space  and  thus  permits  recording  at  fast  tape 
speeds  on  a  small,  low-cost,  low-power  recorder.  A  difficulty 
with  such  a  system,  however,  is  that  the  trigger  is  m  re  likely 
to  be  set  off  by  the  S  wave  or  by  another  later  arrival  than 
by  the  often  smaller  P  wave.  Thus  the  valuable  first 
cycles  of  wave  motion  would  be  lost. 
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The  solution  is  to  build  a  "memory"  into  the  system, 
so  that  the  recorder  can  store  *ta  taken  before  the  threshold 
detector  was  triggered.  One  such  simple  memory  device  is  a 
small,  endless  tape  loop  with  "write"  and  "read"  heads  separa¬ 
ted  by  a  distance  along  the  tape  corresponding  to  the  largest 
probable  S-P  time.  A  second  storage  recorder  is  also  required. 

Such  a  device  tuts  been  used  with  some  success  to  monitor  micro- 
earthquakes  (Omote  et  al .  ,  1955;  Aki  et  al^. ,  1969;  Hunt,  1969). 

A  difficulty  with  such  a  system  in  the  past  has  been  the  severe  in¬ 
crease  in  tape  noise  after  the  thousand  or  more  cycles  necessary 
for  operation  times  of  up  to  a  month.  An  alternative  memory 
device  using  a  single  tape  recorder  and  a  rewind  algorithm  has 
recently  been  described  by  Choudhury  and  Houri  (1973).  Handi¬ 
caps  of  their  system  are  the  regular  periods  when  no  signal 
can  be  recorded  and  a  relatively  short  life  for  the  recorder. 

We  describe  in  this  paper  a  seismometer  system  that  incor¬ 
porates  an  automatic  editor  based  on  a  semiconductor  buffer 
memory  instead  of  a  tape  loop.  The  principles  of  the  editor 
were  briefly  outlined  previously  (Ambuter  and  Solomon,  1972). 

The  semiconductor  memory  has  the  advantage  of  no  moving  parts 
to  degrade  with  repeated  use.  The  recorded  seismic  signal 
is  in  digital  form,  suitable  for  eventual  spectral  analysis. 

The  instrument  package  we  describe,  because  of  our  interest 
in  small  earthquakes,  was  designed  with  a  specific  set  of 
requirements:  (1)  good  response  between  2  and'  30  Hz  frequency; 

(2)  dynamic  range  of  at  least  40  db;  (3)  remote  operation  for 
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periods  of  a  month  or  more;  (4)  low  power;  (5)  low  cost; 

(6)  small  size  and  weight.  Nonetheless,  the  general  charac¬ 
teristics  of  our  system  should  be  applicable  to  a  wide 
range  of  seismic  recording  problems. 

EVENT  DETECTOR 

In  Figure  1  is  shown  a  detailed  block  diagram  of  the 
event  detector  circuitry  used  in  our  system.  Signals  are 
fed  through  a  low-pass  filter  with  a  30  Hz  cut-off  (Figure  2) 
and  are  input  to  a  short-term  averager  (STA)  and,  through  a 
digitally  controlled  switch,  to  a  long-term  averager  (LTA) . 
The  long-  and  short-term  averagers  both  rectify  and  integrate 
the  input  signal.  The  LTA  has  a  time  constant  of  5  minutes, 
which  is  much  longer  than  the  largest  period  of  the  seismic 
signals  we  wish  to  record.  The  LTA  output  is  therefore  re¬ 
presentative  of  the  average  background  level  or  noise.  The 
STA  has  a  time  constant  of  1  second  and  responds  quickly  to 
changes  in  signal  level.  The  output  of  the  LTA  and  a  portion 
of  the  STA  are  continuously  monitored  by  an  event  comparator. 
The  portion  of  the  STA  that  is  fed  to  the  event  comparator 
determines  the  event  threshold.  Put  simply,  an  event  is  de¬ 
tected  when  the  STA  exceeds  the  threshold  level. 

Since  we  are  dealing  with  an  event  recording  system 
with  a  finite  storage  capacity,  one  of  the  prime  considera¬ 
tions  of  the  system  is  to  maximize  the  ratio  of  real  events 
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to  spurious  ones.  Obviously  s  trsde-off  exists  between  the 
threshold  setting  end  the  number  of  events  detected.  Since 
both  the  seismic  signals  of  interest  end  the  noise  may  have 
substantial  energy  in  the  frequency  band  2  to  30  Hr,  restric¬ 
ting  the  event  detection  circuitry  to  look  at  a  portion  of 
this  band  would  probably  not  result  in  improving  the  ratio 
of  detected  real  events  to  .periods  events.  Preliminary 
tests  have  shown  that  for  our  system  a  threshold  of  8  to  10  db 
minimises  the  number  of  spurious  events  while  allowing  all 

real  events  with  signal  to  noise  adequate  for  analysis  to  be 
recorded . 

In  operation,  the  input  to  the  LTA  is  digitally  controlled 
through  a  switch  and  associated  logic.  Under  no-event  condi- 
"ion.,  the  switch  is  in  position  1  (Figure  1)  and  the  LTA 
reflects  the  average  level  of  the  incoming  signals.  When  an 
event  is  detected,  the  switch  is  moved  to  position  2  and  the 
output  of  the  LTA  is  frosen.  The  reason  for  this  is  that  if 
the  LTA  were  allowed  to  respond  to  signal  levels  during  an 
event  in  which  the  signal  amplitude  is  significantly  greater 
than  the  average  background  level,  then  the  LTA  would  charge 
UP  and  cause  premature  termination  of  the  event.  The  LTA 
switch  remains  in  position  2  until  the  STA  falls  below  the 
threshold  level.  The  system  will  continue  to  record  until 
the  STA  level  has  been  below  the  threshold  level  for  a  thirty 

second  time  period.  This  enables  us  to  see  the  trailing  edge 
of  an  event. 


r 
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During  an  event  the  LTA  output  is  fed  into  a  feedback 
amplifier  having  a  gain  such  that  the  slope  of  the  LTA  output 
level  is  slightly  positive.  This  has  minimal  effect  for 
events  of  short  duration,  say  up  to  a  minute.  But  for  events 
which  are  longer  than  typically  expected  from  small  earth¬ 
quake  sources,  the  LTA  will  charge  up  and  eventually  effect 
termination  of  recording.  This  minimizes  the  possibility  of 
exhausting  the  event-recording  storage  capacity  with  data  of 
no  .cientific  interest.  The  LTA  time  constant  was  chosen  to  be 
longer  than  the  longest  signal  period  of  interest,  yet  short 

enough  to  respond  to  noise  build-up,  thereby  minimizing  spurious 
triggering . 


ANALOG  CIRCUITRY 

The  analog  circuitry  (low- noise  amplifiers,  low-pass 
filters,  variable-gain  amplifiers,  sample-and-hold  amplifier, 
variable-gain  logic  control,  and  event  detector)  conditions 
the  incoming  seismic  information  through  amplification  and 
filtering  and  determines  if  an  event  has  been  observed.  The 
digital  circuitry  (A-D  converter,  semi-conductor  memory, 
timing,  and  control  and  format  unit)  samples  the  data  from  the 
output  of  the  analog  circuitry  through  the  A-D  converter  and 
loads  these  data  into  the  semi-conductor  memory  under  the 
direction  of  the  control  and  format  unit.  The  event  detector 
line  is  used  to  turn  on  a  storage  recorder  and  load  the  out¬ 
put  of  the  buffer  memory. 


- n 

-165- 

Seismic  signals  from  the  three-component  seismometer 
(see  Figure  2)  are  first  amplified  using  a  low-noise  integra¬ 
ted-circuit  amplifier ,  and  are  then  passed  through  a  4-pole 
low-pass  active  filter  network  having  a  high-frequency  cutoff 
of  30  Hz.  The  output  from  the  low-pass  filters  is  applied 
to  the  sample-and-hold  circuitry,  and  then  to  the  multiplier, 
which  combines  the  three  channels  into  one  composite  signal. 

One  component  of  data  output  from  the  low-pass  filter,  generally 
the  vertical  component,  goes  to  the  event— detection  circuitry 
previously  described,  and  also  to  the  variable-gain-amplifier 
control  logic. 

The  variable-gain  amflifier  consists  of  an  LTA  and  an 
associated  decoder  blcck.  This  LTA  has  a  time  constant  of 
five  minutes  and  responds  to  the  average  background  or  noise 
level.  The  LTA  output  is  applied  to  the  decoder  block,  whicn 
in  turn  sets  the  gain  of  the  variable-gain  amplifier.  Thus 
the  gain  of  the  system  is  governed  by  the  background  noise 
level.  This  feature  enables  the  seismic  system  to  be  placed 
in  environments  where  the  noise  level  is  unknown  or  subject 
to  variation  and  still  realize  maximum  signal  dynamic  range. 

During  an  event  the  variable-gain  LTA  level  is  frozen  and 
the  system  gain  is  recorded. 

DIGITAL  CONTROL  CIRCUITRY 

The  operation  of  the  remainder  of  the  systen,  governed 
by  the  logic  control  and  format  unit,  may  be  explained  as 


follows : 


1}  The  A-D  converter  upon  command  samples  the  analog 
input  and  generates  an  equivalent  digital  code.  The  system 
being  described  used  an  8  bit  A-D  converter.  ir  operation, 
the  background  noise  is  set  at  approximately  one  half  of  the 
least-significant-bit  by  the  automatic  gain  control,  allowing 
a  signal  dynamic  range  of  42  db. 

2)  The  control  and  format  unit  then  commands  the  A-D 
converter  to  load  its  8-bit  output  words  into  the  semi¬ 
conductor  memory.  The  memory  is  made  up  of'  1000-bit  shift 
registers.  The  input  data  are  shifted  along  until  at  some 
later  time  they  appear  at  the  output  of  the  shift  registers. 
The  capacity  of  the  memory  is  2000  8-bit  words  for  each  data 
channel,  a  total  that  can  easily  be  increased  by  adding  more 
shift  registers.  To  illustrate  the  time  delay,  if  the  sample 
rate  were  200  samples/sec,  then  the  digital  output  would  be 
delayed  10  seconds. 

During  periods  when  the  event  detector  is  off,  the  buffer 
memory  words  just  drop  off  the  end  of  the  shift  register. 

When  an  event  is  detected,  the  storage  recorder  is  turned  on, 
and  the  output  from  the  buffer  memory  is  recorded. 

3)  The  sequence  of  events  that  occurs  when  a  seismic  signal 
is  detected  is  under  the  direction  of  the  control  and  format 
unit.  When  an  event  is  detected,  the  inputs  to  both  long-term 
averagers  are  immediately  inhibited.  This  is  because  the 
long-term  average  is  meant  to  reflect  the 


average  background 


noise,  and  therefore  should  not  be  allowed  to  respond  to 
event  signal  levelr .  Simultaneously,  the  tape  recorder 
motor  is  turned  on  and,  after  a  delay  during  which  the  tape 
motor  comes  up  to  speed,  the  buffer  memory  output  is  re¬ 
corded.  During  the  event,  the  system  gain  and  the  total 
elapsed  time,  expressed  in  terms  of  seconds,  minutes,  hours, 
and  days,  are  recorded  on  the  control  track.  Recording 
continues  until  the  threshold  detector  has  returned  to  its 
presignal  "off"  state  for  a  period  of  30  seconds. 

TIMING  AND  STORAGE  RECORDER 

The  timing  signals  necessary  to  run  this  system  are 
derived  from  a  5  MHz  oscillator.  The  oscillator  output  is 
fed  into  the  time  divider  chain.  Three  sample  speeds  are 
available:  125  Hz,  160  Hz,  and  200  Hz. 

The  storage  recorder  is  a  Sony  800B  four-speed,  portable 
tape  recorder,  modified  by  removing  all  the  audio  circuitry 
and  installing  a  five  track  digital  head.  The  reel  size  is 
five  inches,  allowing  2400  feet  of  recording  tape.  For  a 
200  Hz  sample  rate  and  a  tape  speed  of  1  7/8  inches  per 
second,  a  recording  time  of  four  hours  is  possible.  Propor¬ 
tionately  longer  recording  times  may  be  obtained  by  using 
slower  sampling  rates. 

Three  of  the  recorder's  tracks  are  used  for  storing  data 
from  the  three  input  channels.  The  fourth  track  is  used  for 


recording  of  system  gain  and  associated  timing  marks.  The 
fifth  track  is  for  sync  pulses.  These  sync  pulses  are  used 
to  inform  the  tape  playback  unit  that  there  are  data  on  the 
other  four  tracks  and  also  to  identify  the  relative  bit  posi¬ 
tions  of  the  words  as  they  are  read  from  the  storage  recorder. 
By  using  a  separate  sync  track,  the  system  is  relatively  in¬ 
sensitive  to  tape  speed  variations,  either  in  the  original 
recording  or  in  the  playback  mode.  This  is  a  big  advantage 
over  FM  or  direct  recording  modes,  as  both  of  these  approaches 
are  acutely  sensitive  to  tape  speed  stability  and  impose 
tight  requirements  on  the  storage  recorder. 

The  data  as  recorded  are  not  in  a  format  compatible 
with  large  computers.  The  tape  must  first  be  played  back 
through  a  processor  unit  which  reconverts  the  tape  head  signals 
into  digital  form.  Data  and  control  channels  are  then  applied 
to  the  accumulator  of  a  small  computer  (PDP-7) ,  which  is  pro¬ 
grammed  to  read  the  information  in  the  accumulator  when  the 
sync  channel  produces  a  logic  high.  The  decoded  data  may  then 
be  readily  formatted  and  written  on  an  IBM-compatible  tape. 
Output  from  the  same  processor  unit  may  also  be  fed  into  a 
unit,  containing  D-A  converters,  which  drives  a  multi-channel 
chart  recorder. 

SYSTEM  CHARACTERISTICS 

A  photograph  of  the  system  described  above,  suitable 
for  monitoring  small  earthquakes  on  land,  is  shown  in  Figure 
3.  Included  in  the  figure  are  the  circuit  boards,  storage 


recorder,  and  geophones  that  constitute  the  functional  part 
of  the  system.  All  instrumentation  aside  fr  m  sensors  and 
batteries  can  fit  into  a  suitcase  of  dimensions  9"  x  10"  x  12" 
for  land  use  into  a  9  1/2"  diameter  I.D.  cylinder  for  ocean- 
bottom  experiments.  The  weight  of  the  system,  again  excluding 
sensors  and  batteries,  is  about  twenty  pounds. 

System  power  requirements  are  dependent  on  the  stability 
of  the  timing  oscillator.  When  using  relative  arrival  times 
of  P  and  S  waves  at  array  stations  to  locate  small  earthquakes, 
a  relative  time  accuracy  of  at  least  .01  sec  is  often  desirable. 
Por  land  stations  in  areas  where  WWV  radio  can  be  received,  a 

7 

timing  stability  of  1  part  in  10  per  day  is  sufficient.  The 
power  drain  for  a  system  such  as  we  describe  with  this  timing 
stability  averages  about  three-quarters  of  a  watt.  For  1-month 
recording  times  at  stations  remote  from  an  external  timing 
reference,  such  as  on  the  ocean  bottom,  an  oven-controlled 

9 

oscillator  with  timing  stability  of  1  part  in  10  per  day  is 
necessary.  The  system  power  drain  is  then  about  1  1/2  watts. 

The  voltages  required,  12  volts  and  ±5  volts,  may  be 
supplied  by  many  types  of  batteries.  We  have  found  both  lead- 
acid  batteries  and  magnesium  batteries  to  be  useful.  Lead- 
acid  batteries  are  rechargeable,  but  are  comparatively  bulky. 
Magnesium  batteries  are  more  expensive,  but  are  much  lighter. 
Thirty  days  of  system  operation  require  about  $120  worth  of 
either  type  of  batteries. 
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The  total  component  cost  of  the  system  shown  in  Figure  3, 
including  geophones  (4.5  Hz  frequency),  is  approximately 
$2000  with  the  oven-controlled  oscillator,  and  is  $200  1  »ss 
with  the  less  demanding  timing  requirements.  About  two  man- 

weeks  of  construction  time  are  necessary  to  assemble  the 
complete  system. 


SYSTEM  TESTS 


The  event-recording  circuitry  was  tested  using  data 
from  the  Agassiz  seismic  observatory  at  Harvard,  Massachusetts. 
The  system  sampling  rate  was  160  samples/sec,  the  buffer 
memory  delay  was  12.5  seconds,  and  the  maximum  signal  fre¬ 
quency  was  20  Hz.  The  signal  from  a  short-period  vertical 
seismometer  was  processed  and  output  to  a  two-channel  heli- 
corder.  The  first  channel  of  the  helicorder  had  as  its  input 
the  continuous  output  of  the  system's  buffer  memory.  Input 
to  the  second  helicorder  channel  was  derived  from  the  buffer 
memory  output  gated  with  the  output  of  the  event  detector 
logic  (which  turns  the  storage  recorder  on  and  off) .  in  such 
an  operation,  the  first  helicorder  trace  is  a  continuous 
representation  of  the  incoming  ihort-period  data.  The  second 
♦ra  e  is  a  straight  line  unless  the  signal  threshold  is  ex¬ 
ceeded,  in  which  case  the*  two  helicorder  traces  are  identical. 

The  threshold  level  was  initially  set  3  ob  and  gradu¬ 
ally  increased  to  a  maximum  value  of  14  db. 


The  number  of 
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event  triggerings  per  unit  time  depended  strongly  on  the 
threshold  level.  At  a  detection  threshold  of  6  db,  we 
recorded  about  60  to  80  "events"  per  day,  though  the  great 
majority  were  clearly  noise  rather  than  small  earthquakes 
or  explosions.  At  a  detection  threshold  of  12  db,  we  recorded 
real  seismic  events  only,  less  than  1  per  day,  though  a  few 
small  signals  from  possible  sources  of  interest  were  missed. 

An  optimum  threshold  level  for  this  site  lies  probably  be¬ 
tween  8  and  10  db.  At  this  seating,  the  number  of  spurious 
triggers  was  minimized  without  missing  any  clearly  identifiable 
events.  At  a  different  site,  for  instance  one  less  noisy, 
the  optimum  threshold  setting  might  be  different. 

Portions  of  helicorder  recordings  shown  in  Figure  4 
illustrate  typical  detected  events.  The  operation  of  the 
buffer  memory  delay  is  clearly  illustrated  in  the  figure; 
both  noise  prior  to  each  event  and  the  unambiguous  first 
arrivals  have  been  recorded.  The  system  continued  to  record 
until  the  STA  stayed  below  the  threshold  level  for  a  period 
of  30  seconds. 


CONCLUSIONS 

We  have  described  a  seismic  recording  system  that  is 
ideal  for  monitoring  small  earthquakes  in  remote  areas. 

The  system  is  small,  consumes  little  power,  and  is  low  in 
cost.  This  is  in  large  part  due  to  the  latest  advances  in 


integrated-circuit  technology.  For  example,  the  2000-word 
8-bit  memories  in  this  system  cost  less  than  $100  and  occupy 
an  area  of  about  2  square  inches  on  a  circuit  board.  The 
dynamic  range  is  large  and  may  be  readily  extended  by  using 
a  10-  or  12-biu  A-D  converter. 

The  system  is  insensitive  to  tape-speed  variations, 
which  are  critical  to  other  recording  techniques.  In  play- 
back  mode,  all  that  Is  required  is  to  distinguish  between  a 
logic  high  and  a  logic  low.  Therefore  the  specifications  of 
the  storage  recorder  are  not  particularly  stringent  and  a 
small,  inexpensive  reoorder  may  be  used. 

This  system  may  be  readily  modified  to  meet  different 
applications.  Sampling  rate  may  be  varied  by  a  factor  of  up 
to  10  simply  by  changing  one  wire.  Additional  data  channels 
may  be  added  at  a  cost  of  $150  per  channel.  In  applications 
where  long-period  data  ar*.  of  interest,  the  operation  time 
could  easily  be  extended  to  a  year.  It  should  be  noted  that 
the  sensor  need  not  be  a  geophonc?  thus  the  above  system  could 
be  applied  to  other  fields  with  minimal  functional  changes. 

By  design  this  system  is  self-editing,  so  that  data 
analysis  does  not  need  to  be  preceded  by  a  tedious  search  of 
voluminous  tapes  for  useful  data.  The  signals  are  recorded 
in  digital  form,  suitable  for  eventual  spectral  analysis  on  a 
digital  computer. 

The  details  of  the  system  described  were  designed  for  the 
specific  aim  of  monitoring  small  earthquakes  on  the  ocean- 
bottom  or  in  remote  land  areas  of  relatively  low  seismicity. 
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Nonetheless,  the  requirements  met  by  such  a  system  are 
shared  to  some  extent  by  such  diverse  experiments  as  seismic 
recording  by  unmanned  planetary  landers  (Anderson  et  al .  * 
1972}  ,  unmanned  seismic  observatories  of  a  long-term  or 
permanent  nature,  and  even  routine  data  processing  in  real¬ 
time  of  a  conventional  seismometer  network.  Instruments  of 
similar  design  should  see  widespread  implementation  by 
seismologists  in  the  next  few  years. 


-174- 


ACKNOWLEDGEMENTS 

We  thank  Keiiti  Aki,  Paul  Reasenberg  and  Paul  Mattaboni 
for  many  helpful  suggestions  during  the  design  and  develop¬ 
ment  stages  of  this  work. 

'iMs  research  was  partially  sponsored  by  the  Advanced 
Research  Projects  Agency  and  monitored  by  the  Air  pcrce 
Office  of  Scientific  Research  under  contract  F  44620-71-C-0049 . 


REFERENCES 


Aki ,  K. ,  M.  Hori ,  and  H.  Matumoto  (1969) .  Microftf tershocks 
observed  at  a  temporary  array  station  on  the  Kenai 
Peninsula  from  May  19  to  June  7,  1P64,  Coast  and  Geodetic 
Survey  Publication  10-3,  131-156. 

Ambuter ,  B.P.,  and  S.C.  Solomon  (197P)  .  An  automatic  seismic 
editing  system  (abstract) ,  EOS  Trang.  Amer.  Geophys .  Un. , 
53,  105.1. 

Anderson,  D.L.,  R.L.  Kovach,  G.  Latham,  F.  Press,  M.N .  Toksflz, 
and  G.  Sutton  (1972).  Seismic  investigations:  The 
viking  Mars  Lander,  Icarus  16 ,  205-216. 

Choudhury ,  M.A.,  and  A.  Hour!  (1973).  A  low-cost  observatory 
tape  recorder,  Bull.  Seism.  Soc.  Am.  63,  877-884. 

Dibble,  R.R.  (1964)  .  A  portable  slow  motion  magnetic  tape 

recorder  for  geophysical  purposes,  N.Z.J.  Gcol .  Geophys. 
7,  445-465. 

Green,  R.  (1973).  A  portable  multi-channel  seismic  recorder 
and  a  data  processing  system.  Bull.  Seism.  Soc.  Am..  63 , 
423-431. 

Hunt,  D.P .  (1969)  .  A  selected  events  recorder  for  seismic 
applications,  Mackay  School  of  Mines  Technical  Report, 
University  of  Nevada,  Reno,  141  pp. 


i 


-176- 


Muirhead,  K.J.,  and  D.W.  Simpson  (1972).  A  three-quarter 
watt  seismic  station,  Bull.  Seism.  Soc.  Am.  62, 

985-990. 

Omote,  S.,  S.  Miyamura,  and  Y.  Yamazaki  (1955).  Triggered 
magnetic  tape  recorder  for  routine  seismic  observation, 
Bull.  Earthquake  Res.  Inst..  Univ.  Tokyo  33.  397-40? 

Sacks,  I.S.  (1966) .  A  broad-band  large  dynamic  range  seismo¬ 
graph,  in  The  Earth  Beneath  the  Continents,  ed.  by 
J.S.  Steinhart  and  T.J.  Smith.,  Amer.  Geophys.  Un.  Mon. 

10,  543-553. 


-177- 


PI  GURE  CAPTIONS 

Figure  1.  Block  diagram  of  event  dutector  circuitry. 

Figure  2.  Block  diagram  of  complete  seismic  recording  system. 

Figure  3.  Photograph  of  system  suitable  for  remote  operation 
o.i  land,  disassembled  for  clarity. 

Fic/ure  4.  Examples  of  signals  which  triggered  event-detection 
circuitry.  Source  are  probably  quarry  blasts.  Times 
are  in  G.M.T.  Tic  marks  are  generally  spaced  at  1  minute 
intervals . 
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